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ABSTRACT 

Breather  solitons  in  a  one-dimensional  lattice  of  coupled  nonlinear  oscillators  are 
numerically  investigated.  These  are  localized  nonpropagating  steady  states  that  exist  at 
frequencies  either  below  the  linear  cutoff  frequency  (corresponding  to  the  extended  mode  in 
which  all  the  oscillators  are  in-phase)  or  above  the  upper  linear  cutoff  frequency  (corresponding 
to  the  extended  mode  in  which  each  oscillator  is  180°  out-of-phase  with  its  immediate  neighbors). 
The  lattice  is  damped  and  parametrically  driven.  A  nonlinear  Schrodinger  theory,  which  assumes 
a  modulational  amplitude  that  is  weakly  nonlinear  and  slowly  varying  in  space,  is  compared  to 
the  numerical  data.  The  error  is  roughly  5%  at  low  amplitudes  and  20%  at  high  amplitudes. 
The  regions  in  the  drive  parameter  plane  (amplitude  vs  frequency)  where  the  breathers  exist  are 
numerically  determined  and  compared  to  theory.  A  substantial  discrepancy  occurs  at  lower  drive 
amplitudes  where  the  theory  predicts  that  the  lower  cutoff  breather  should  exist,  but  where  an 
instability  is  observed.  Also  in  contrast  to  the  theory,  the  region  of  the  upper  cutoff  breather  has 
relatively  large  areas  in  which  quasiperiodicity  occurs  or  the  motion  decays  to  rest. 
Quasiperiodicity  is  also  observed  in  the  lower  cutoff  breather.  Finally,  instead  of  a  global 
parametric  drive,  an  end  drive  is  investigated.  It  is  found  that,  for  drive  frequencies  outside  the 
linear  propagation  band,  there  is  an  amplitude  threshold  for  the  periodic  ejection  or  "shedding" 
of  propagating  breather  solitons.  The  quasiperiodicity  that  occurs  for  a  global  parametric  drive 
may  be  a  consequence  of  soliton  shedding. 
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I.  INTRODUCTION 

A.       BACKGROUND 

A  soliton  is  a  localized  wave  whose  shape  or  envelope  is  constant  as  a  result  of 
nonlinearities  and  dispersion.  Furthermore,  it  collides  elastically  with  other  waves.  Solitons 
were  first  observed  as  canal  surface  waves  by  John  Scott  Russell  in  1834  (Dodd  et  al.  1978). 
Other  systems  are  now  known  to  support  solitons  of  various  types.  Research  in  fiber  optic 
solitons  is  currently  very  active  because  of  their  potential  use  in  communications  and  optical 
switching  (Mollenauer  et  al.  1991). 

In  this  thesis,  a  one-dimensional  lattice  of  coupled  nonlinear  oscillators  is  numerically 
investigated.  This  system  is  known  to  possess  kink  solitons  (Galvin  1990,  Denardo  et  al.  1991). 
It  will  be  shown  that  breather  solitons  can  also  exist.  These  are  nonpropagating  steady  states  in 
which  most  of  the  oscillators  are  approximately  at  rest  while  those  in  a  localized  region  have 
substantial  amplitude.  Such  a  state  is  clearly  not  a  linear  mode  of  the  system.  Breathers  are  self- 
trapped  states  that  occur  at  frequencies  outside  the  linear  propagation  band.  Two  types  exist: 
lower  cutoff  breathers,  in  which  the  oscillators  are  in-phase,  and  upper  cutoff  breathers,  in  which 
the  oscillators  are  180°  out-of-phase.  Lower  and  upper  cutoff  solitons  can  be  considered  as 
finite-amplitude  modulations  of  the  uniform  lower  and  upper  cutoff  modes,  respectively. 

In  the  absence  of  drive  and  dissipation,  it  is  known  that  weakly  nonlinear  breathers 
described  by  the  </>"  theory  are  unstable,  although  the  decay  rate  is  extremely  slow  (Segur  and 
Kruskal  1987).  It  is  also  known  that  instabilities  can  occur  as  a  result  of  discreteness  in  systems 
that  are  undriven  and  undamped  (Goedde  1990).  In  this  thesis,  dissipation  and  global  drive  are 
employed,  which  lead  to  stable  breathers. 
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The  simplicity  of  the  model  equation  suggests  that  breathers  can  occur  in  a  variety  of 
systems.  Indeed,  the  first  observation  of  breathers  was  as  cross  surface  waves  on  a  long  channel 
of  liquid  (Wu  et  al.  1984).  Lower  cutoff  breathers  have  been  observed  in  a  pendulum  lattice 
(Denardo  1990).  Recently,  upper  cutoff  breathers  have  been  observed  in  a  magnetically  coupled 
pendulum  lattice  and  confirmed  numerically  with  an  equation  that  models  this  system  (Atchley 
1991). 

It  is  also  shown  that,  for  a  local  pure-frequency  drive  confined  to  one  end  site  of  a  lattice, 
it  is  possible  to  shed  propagating  solitons.  The  shedding  occurs  at  a  frequency  that  is 
quasiperiodic  relative  to  the  drive  frequency.  This  phenomenon  has  previously  been  observed 
in  only  one  system,  surface  waves  just  below  the  second  cutoff  mode  in  a  large  tank  (Kit  et  al. 
1987,  Shemer  1990).  The  existence  of  the  soliton  shedding  in  a  simple  lattice  suggests  that  it  is 
a  general  phenomenon  which  can  occur  in  a  variety  of  systems,  and  that  a  simple  fundamental 
explanation  should  exist,  although  no  such  explanation  is  known  at  present. 

B.       EQUATION  OF  MOTION 

To  motivate  the  equation  of  motion  that  will  be  investigated,  a  model  system  is  considered. 
The  system  is  a  uniform  lattice  of  linearly  coupled  simple  pendulums  that  oscillate  transverse  to 
the  lattice.   The  torque  on  the  nth  pendulum  is: 

rn  -  /x   (  0n+i   -  20n    +  0n_!   )    -  mgl  sin(0n),  J.B.I 

where  \i  is  the  torsional  constant  that  characterizes  the  coupling,  m  is  the  pendulum  mass,  and 
1  is  the  pendulum  length.   For  weakly  nonlinear  motion,  the  approximation 

sin(0)     =  0   -  -  03  J.B.2 

v    '  6 


sin(0)    =  6   -  1  03  I.B.2 

v   '  6 

can  be  made.  Linear  damping  is  included  by  adding  to  LB.  1  a  torque  proportional  to  the  velocity 
d'D  (with  a  negative  coefficient),  where  the  prime  denotes  differentiation  with  respect  to  time. 
To  sustain  localized  structures  in  the  lattice,  a  global  parametric  drive  that  results  from  vertically 
oscillating  the  lattice,  is  utilized.  By  the  Equivalence  Principle,  in  the  reference  frame  of  the 
lattice  the  effective  acceleration  due  to  gravity  is: 

9eU  -  g  +  a  cos(2cot)  ,  I.B.2 

where  a  is  the  acceleration  amplitude  of  the  drive  and  2w  is  the  drive  frequency. 

Combining  the  above  effects,  and  using  Newton's  Second  Law,  the  equation  of  motion  of 
for  the  lattice  can  be  expressed  as: 

8'i-c2(8n+1  -  26n  +  Qn_x)+W  n+[u02+2y\cos(2ut)]Qn  -  a6n3,    I.B.4 

where  c2  is  a  measure  of  the  coupling  between  lattice  sites  and  co0  is  the  cutoff  frequency  of  a 
pendulum.  For  simplicity,  the  term  proportional  to  the  product  of  the  drive  and  nonlinearity  has 
been  dropped.  It  can  be  shown  that  this  has  no  essential  effect  upon  the  localized  structures.  The 
nonlinear  coefficient  a  equals  a>02/6  if  (IB. 4)  is  to  approximate  a  pendulum  lattice.  The  model 
is  generalized  by  allowing  the  nonlinear  coefficient  to  be  positive  or  negative.  The  lattice  is 
driven  parametrically  with  an  amplitude  of  2t\.  The  fundamental  response  frequency  u  is  half  the 
drive  frequency.  The  lower  cutoff  mode  of  the  lattice  is  characterized  by  uniform  motion  in  the 
lattice  (Fig.  I.B.I),  and  has  linear  frequency  co0.  The  upper  cutoff  mode  has  a  180°  phase 
difference  between  adjacent  lattice  sites  (Fig.  I.B.2),  and  has  linear  frequency  (oo02  +  4c2)"2.  For 


a  softening  (a  >  0)  lattice,  the  upper  cutoff  mode  is  stable.  The  lower  cutoff  mode  is  subject 
to  the  Benjamin-Feir  instability  (Denardo  1990),  and  the  motion  evolves  into  one  or  more 
breather  solitons  (Fig.  I.B.3).  For  a  hardening  (a  <  0)  lattice,  the  roles  of  the  cutoff  modes  are 
reversed  and  an  upper  cutoff  breather  evolves  (Fig.  I.B.4). 


Figure  I.B.I   Lower  Cutoff  Mode 


Figure  I.B.2   Upper  Cutoff  Mode 
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Figure  I.B.3   Lower  Cutoff  Breather 


! 

• 

. 

A                    A                  A                    A 

Figure  I.B.4   Upper  Cutoff  Breather 


II.  THEORY 

A.   LOWER  CUTOFF  BREATHER 

For  modulations  of  the  lower  cutoff  mode,  the  displacement  (0J  in  I.B.4  can  be  written 
approximately  as: 

d n-A(n,  t)  eiv>t+  complex  conjugate,  II. A. 1 

where  A(\,t)  is  a  complex  differentiate  function  of  its  arguments.  The  lattice  spacing  is  assumed 
to  be  unity.  A  displacement  (0J  containing  a  third  harmonic  term  produced  approximately  the 
same  results  as  II. A.  1  and  was  consequently  discarded.  Substituting  II. A.  1  into  I.B.4  and 
requiring  that  A  be  slowly  varying  and  weakly  nonlinear,  one  obtains  the  modulational  equation 

2ioA  -  c2  A^  +    (0)o2  -  O)2  +  iu>P)A  +  r|  A*  -  3a  I  A  I2  A,  II. A. 2 

which  is  a  nonlinear  Schrodinger  equation.  If  the  modulation  is  assumed  to  be  nonpropagating, 
a  stationary  solution 

A   (x,  t)    -  A   (x)    eia  XT. A. 3 

exists  where  A  and  6  are  real.   Substituting  II. A. 3  into  II. A. 2  gives: 

-  c2  A"  +    (co02  -  G>2)  A  -  3a  A3  +  ia>PA  -  -TiAe'2i8.  IJ.A.4 

Equating  the  imaginary  parts  of  both  sides  yields  an  expression  for  5: 

Sin(26)    -  -^-.  JJ.A.5 


The  real  part  of  II. A. 4  yields: 


c2  A"  -   \i  A  +  3   a  A3  -  0, 
where 


II. A. 6 


[X  -  a>  02  -  (D2  +  v^l2  -  <*>2P: 


There  are  actually  two  possibilities  because  of  the  radical.    However,  for  a  softening  lattice  the 
negative  branch  leads  to  an  unstable  solution  (Denardo  1990). 

Equation  II. A. 6  can  be  integrated  to  yield  an  elliptic  function.    A  localized  solution  to 
II. A. 6  is  of  the  form  A  =  a  sech(bx),  provided 


2c2 
NlTo" 


and     b 


\ 


JL 


Then 


A  - 


3a 


sech 


\ 


-%    U  -  x0) 

c2 


II. A. 1 


Substituting  the  solution  for  A  into  II.A.l,  the  expression  for  the  displacement  is: 


6 


N 


Hi 

3a 


sech 


(x 


*o> 


cos  (a)  t  +  b) 


II. A. 8 


This  solution  is  valid  for  /*  >  0  and  a  >  0  (softening). 

The  physical  reasoning  for  the  existence  of  the  breather  is  derived  from  the  curvature  of 
the  modulation  envelope.  For  the  softening  lower  cutoff  breather,  negative  curvature  in  the  body 
of  the  envelope  produces  a  restoring  force  (Fig.  II.A.l).  The  frequency  in  the  body  is  below 
cutoff  due  to  the  nonlinearity  which  decreases  the  frequency  more  than  the  curvature  increases 
it.  In  the  tail  of  the  envelope  the  curvature  is  positive,  thus  an  anti-restoring  force.  The 
frequency  is  below  cutoff  and  the  breather  is  evanescent  in  the  tail.    At  the  inflection  points,  the 


nonlinearity  is  the  sole  reason  the  frequency  is  below  cutoff.  The  breather  is  thus  a  self-trapped 
state  due  to  the  combination  of  the  nonlinearity  and  curvature  of  the  modulation.  As  a  result, 
steady  state  motion  is  possible.   This  motion  has  been  observed  in  the  simulated  lattice. 

Concluding  that  steady  state  motion  is  possible,  one  should  be  able  to  determine  the  values 
of  7}  and  o)  where  a  stable  solution  exists.   Requiring  /x  from  II. A. 6  to  be  positive, 

x\2  >    (   0)o2  -  O)2   )2  +  g>2  P2.  II. A. 9 

If  u~oi0,  this  yields  the  Matthieu  hyperbola.  The  sech  solution  (Eqn.  II. A. 8)  is  unstable  inside 
the  hyperbola  because  of  growth  in  the  wings.  There  are  two  conditions  for  the  solution  to  exist: 
1}  >  co/3  and  fi  >  0.  The  first  condition  is  i\  >  ujt  if  co«co0.  The  second  condition  is  related 
to  II. A. 9,  which  can  be  rewritten  as 

si  r\2  -  G)2  P2     >  lto02  -  o)2l.  JJ.A.10 

If  w  >  co0,  then  y.  >  0.  However,  if  w  <  co0  the  requirement  for  y.  >  0  is  77  >  co/8.  The  sech 
solution  decays  to  rest  below  the  boundary  77  >  co/3. 
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Figure  I I.A.I   Lower  Cutoff  Breather  Envelope  Curvatures 


B.       UPPER  CUTOFF  BREATHER 

The  displacement  (0J  for  modulations  of  the  upper  cutoff  mode  can  be  written  as 

Qn-   (-1)  n  A(n,  t)  eiwt+  complex  conjugate  II.B.l 

Where  A(n,i)  is  a  complex  differentiable  function  of  its  arguments.  Again,  the  higher  order  terms 
have  been  neglected,  and  the  lattice  spacing  is  assumed  to  be  unity.  The  equation  that  describes 
the  weakly  nonlinear  complex  amplitude  is  again  a  nonlinear  Schrodinger  equation: 

2i(jiAt  +   c2  A^  +    (to!2  -  a)2  +   ioP)A  +  T|  A  -   3a  I  A  I2  A.  II.  B.  2 

Equation  II. B. 2  is  identical  to  II. A. 2  with  the  substitutions  w02-»w,2  and  c2-*-c2,  where  a*  =  w02 
+  4c2,  which  is  the  square  of  the  linear  frequency  of  the  upper  cutoff  mode.  Assuming  a 
nonpropagating  modulation,  a  stationary  solution  for  A  is: 

A    (x,  t)    -  A    (x)    eih ,  II.  B.  3 

where  A(x)  and  8  are  real.   Substituting  II. B. 3  into  II. B. 2  gives: 

c2  A"  +    (u>2  -  w2)A  -  3a  A3  +  ia>PA  -  -x]Ae'2it> .  II.  B.  4 

Setting  the  imaginary  parts  of  both  sides  equal  and  solving  for  6,  the  phase  relation  is  again: 

sin(26)    -  -^£.  II.  B.  5 

n 

The  real  part  of  II. B. 4  yields: 

-c2  A"  +  v  A  +  3   a  A3  -  0  II.  B.  6 

Where  v  -  ca2  -  w^  +  \]r\2  -  to2P2 
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which  is  similar  to  II. A. 6.    Again,  the  selected  sign  of  the  radical  leads  to  a  stable  solution.   A 
trial  solution  of  the  form  A  =  a  sech(bx)  used  in  II. B. 4  yields: 


A  - 


2v 


\~^ 


sech 


\ 


c2 


II.B.l 


Substituting  the  solution  for  A  into  II.B.l,  the  displacement  (0J  for  the  upper  cutoff  mode  is: 


en  -  (-D  n  2 


2v 


\^ 


sech 


\ 


—    (x  -  x0) 

c2 


cos  (cat  +  6)  ,     II. B.  8 


which  is  valid  for  »>>0  and  a<0  (hardening). 

The  physical  mechanism  for  the  existence  of  the  upper  cutoff  breather  is  essentially 
identical  to  that  of  the  lower  cutoff  breather.  The  negative  curvature  in  the  body  of  the  hardening 
upper  cutoff  breather  now  produces  an  anti-restoring  force  (Fig.  II.B.l)  (Denardo  1990), 
however,  the  nonlinearity  is  stronger  than  the  curvature  and  the  frequency  in  the  body  is  above 
the  linear  upper  cutoff  frequency.  The  positive  curvature  in  the  tails  of  the  modulation  creates 
a  restoring  force  and  the  frequency  is  again  above  the  cutoff  frequency.  Similar  to  before,  the 
frequency  is  above  the  cutoff  frequency  at  the  inflection  points  as  a  sole  result  of  the  nonlinearity. 
Thus,  the  upper  cutoff  breather  is  also  a  self-trapped  state  due  to  the  combination  of  curvature 
and  nonlinearity. 

As  for  the  lower  cutoff  case,  a  region  of  stable,  steady  state  motion  should  also  exist  for 
the  upper  cutoff  breather.   Constraining  v  from  II. B. 6  to  be  positive, 


il2  >    (    O),2  -   w2   )2   +  co2   P2. 


II.  B.  9 


11 


Again,  if  co  =  co0  this  yields  the  Matthieu  hyperbola.  The  wings  of  the  sech  solution  are  unstable 
inside  the  hyperbola.  There  are  two  conditions  for  the  solution  exist:  r\  >  co/3  and  v  >  0.  The 
first  condition  is  77  >  co<jS  if  co  =  co0.  The  second  condition  is  related  to  II. A. 20  which  can  be 
rewritten  as 

7  Tl2  "  w2  P2      >  l^i2  "  w2l-  JJ.B.10 

If  co  >  co,,  then  n  >  0.  However,  if  co  <  w,  the  requirement  for  ^t  >  0  is  17  >  co/S.  The  sech 
solution  decays  to  rest  below  the  boundary  77  >  co/3. 
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Figure  II.B.l   Upper  Cutoff  Breather  Envelope  Curvatures 
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C.       END  DRIVEN  LATTICES 

Nonlinear  lattices  driven  at  one  end  and  at  a  frequency  outside  the  linear  propagation  band 
have  been  observed  to  periodically  eject  solitons  (Sect. III. D).  This  is  in  strong  contrast  to  the 
linear  theory  which  predicts  only  a  steady  state  evanescent  wave,  as  will  now  be  shown.  The 
equation  of  motion  for  a  lattice  of  linearly  coupled  nonlinear  oscillators  (I.B.4)  can  be  modified 
to  replace  the  parametric  drive  with  a  direct  drive.  Further  restricting  the  drive  to  the  end  site 
alone  and  removing  the  periodic  boundary  condition  yields: 

23-0:       e/0/-c2(61  -  60  )+p6/0  +  (0o2  60+Ticos(o)t)    -  cc603  II.  C.  2 

n>o:      e'i-c2^  -  en)+pev«e2  °*  -  «°„3- 

This  produces  an  end  driven  lattice.  The  lattice  is  still  linearly  coupled  and  nonlinear,  but  is  now 
more  suitable  for  the  observation  of  propagating  solitons  which  have  a  spatially  varying  phase 
(Denardo  1990).  Focusing  on  the  linear  motion  in  the  continuum  limit,  equations  II. C. 2  have 
the  form: 

x  >   0    :  ytt  -  c2  Yxx  +   P  yt  +  fa)  02  y  -  0  JJ.  C.  3 

c2 
x  -  0    :  ytt yx  +  <j)  02  y  +  r\  cos  (cot)  . 

G. 

Consider  a  solution  of  the  form: 

y  -  A  ei*x_iut"ffx  +  complex  conjugate,  II.  C. 4 

with  a  complex  wavenumber  K  =  k  +  ia.   Substituting  II. C. 4  into  II. C. 3  yields: 

-o)2  -  c2   (ik  -  a)2  -  iop  +  o>  2  -  0 .  II.  C.  5 
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Separating  the  real  and  imaginary  parts  of  II. C. 5  gives: 


G)    2    -    G)2 

Re:  a2  -  k2 °- 


J/n: 


2C2 


The  case  of  interest  is  o>  <  co0.   The  solution  for  /3=0  is: 


k-0        and       a 


N 


G)  02    -    O)2 


II.  C.  6 


Physically,  the  wave  is  purely  exponential  in  space;  the  sign  of  the  displacement  is  the  same  for 
every  lattice  site.  If  /S^O  then  k  can  be  eliminated  using  the  real  and  imaginary  parts  of  II. C. 5 
to  obtain: 


O)    2    -    G)2 
ff4     _    2     J^O ^_     a2     _ 

2c2 


'    G)P 

,  2c2 


II.  C.  7 


Solving  for  a2  and  choosing  the  positive  root ,  since  a  was  assumed  to  be  real,  gives: 


tt2 1-    [tt  02    -    G)2    +    J(G)Q2    -    G>2)2    +     (G)p)2l 

2c2     L  ¥  J 


JJ.C.8 


The  solution  for  k  is  then: 


Jc-  J^A 

2c2   a 


JJ.C.9 


If  the  solutions  for  u>2  and  k  are  specialized  to  weak  damping,  w|8  <  co02-co2  then: 


a  - 


N 


G)  Q2    -    G)2 


2c    ^ 


JJ.C.10 


2    -    G)2 
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The  clamping  gives  rise  to  a  spatial  phase  (k  ^  0).  Specializing  to  strong  damping,  co/S  >  <o02-w2 
or  near-propagation,  yields: 


a  -  k 


\    2, 


a>P 


II.  C.  11 


This  is  the  largest  that  k/a  can  be  (unity).   In  general,  0<k/a<  1  (Fig.  II.C.l). 
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Figure  II.C.l   Evanescent  Wave  For  k/a=l 
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III.   NUMERICAL  SIMULATION 

A.       IMPLEMENTATION 

The  motion  of  the  lattice  of  linearly  coupled  nonlinear  oscillators  was  solved  with  the 
Euler-Cromer  method.  This  implementation  utilized  acceleration  as  the  basic  numerical  iteration. 
The  position  calculation  of  each  lattice  site  after  a  given  time  step  was  divided  into  three  sections: 

1 .  The  calculation  of  acceleration  for  each  lattice  site  due  to  the  position  of  it  and  the 
two  adjacent  sites. 

2.  The  calculation  of  the  new  velocity  for  each  site  by  multiplying  the  acceleration  by  the  time 
step,  and  adding  the  old  velocity. 

3.  The  calculation  of  the  new  position  for  each  lattice  site  using  the  new  velocity  for  each  site. 
Periodic  boundary  conditions  were  imposed,  thus  making  the  lattice  a  ring  lattice  with  any 
curvature  ignored.  The  Euler-Cromer  method  converged,  with  an  increasing  number  of  time 
steps,  to  a  "true"  lattice  amplitude  (Fig.  III. A.  1).  A  time  step  that  gave  an  amplitude  within  .  1  % 
of  the  "true"  response  amplitude  allowed  the  program  to  run  at  a  fairly  high  rate  of  speed  with 
minimal  loss  of  accuracy. 

The  program  is  highly  interactive  and  gives  a  "real  time"  picture  of  the  lattice  motion.  The 
numerical  simulation  imposes  some  interactive  restrictions.  The  elapsed  time  in  the  model  frame 
must  be  reset  periodically  (preferably  every  period  of  the  drive)  to  avoid  strong  transients  when 
changing  the  drive  frequency.  The  problem  stems  from  the  parametric  drive  term  (cos(2cot)) 
since  the  time  increment  is  based  on  the  response  period.  Thus  a  change  in  frequency  in  the 
middle  of  the  cycle  produces  a  discontinuity  in  the  size  of  the  time  increment.  The  solution  is 
to  employ  an  integer  counter  to  reset  the  elapsed  time  after  one  period  of  the  drive  and  to  require 
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any  frequency  changes  to  occur  when  the  time  is  reset.  The  startup  of  a  modulation  envelope  or 
"profile"  can  be  accomplished  in  two  ways:  restarting  a  saved  profile,  or  starting  a  profile  created 
from  theory.  A  profile  can  be  saved  while  the  program  is  running,  thus  eliminating  the  need  to 
duplicate  steps  required  to  reach  a  particular  set  of  parameters.  All  necessary  parameters  such 
as  drive  amplitude  and  frequency  in  addition  to  the  position  and  velocity  for  each  lattice  site  are 
written  to  a  file  of  the  researcher's  choice.  The  modulation  must  meet  certain  stability  criterion 
to  be  recorded  and  is  "captured"  at  the  upper  turning  point  (response)  of  a  selected  site.  A 
previously  saved  profile  can  be  recalled  using  only  the  filename.  Since  the  response  and  drive 
have  a  phase  difference,  the  program  restarts  the  profile  with  a  time  phase  correction  to  minimize 
startup  transients.  A  startup  file  could  also  be  generated  using  one  of  the  programs  that  calculate 
a  profile  from  theory  using  keyboard  input  of  drive,  coupling  and  other  parameters.  The 
equations  of  motion  are  in  a  subroutine,  allowing  them  to  be  altered  or  replaced  quickly.  A 
different  method  of  calculating  the  lattice  site  positions  would  probably  not  cause  any  difficulty 
as  long  as  the  new  coordinates  were  available  at  the  same  place  in  the  program  as  before.  A 
complete  list  and  more  detailed  description  of  the  program  options  as  well  as  a  program  code 
listing  are  included  in  the  Appendix. 
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Figure  III.A.l   Response  Accuracy  Plot 


20 


B.       LOWER  CUTOFF  BREATHERS 

In  order  to  investigate  the  stability  of  lower  cutoff  breathers,  a  drive  plane,  consisting  of 
drive  amplitudes  and  frequencies  where  the  particular  structure  exists,  was  obtained  by  numerical 
solution  of  the  equation  of  motion.  The  structure  chosen  for  the  lower  cutoff  drive  plane  was  a 
symmetric  breather  (Fig.  III.B.l)  in  which  the  maximum  amplitude  is  centered  on  a  lattice  site 
(on-site).  The  nonlinear  parameter  (a)  was  1  (lower  cutoff),  the  coupling  (c2  or  7)  was  .1 
(weak),  dissipation  (/3)  was  .03  (small)  and  the  lattice  consisted  of  50  sites  long.  The  points  in 
the  drive  plane  (Fig.  III.B.2)  were  almost  exclusively  obtained  by  incrementing  the  drive 
amplitude  by  .1%  (or  less)  and  recording  the  amplitude  at  which  transitions  occurred. 
Additionally,  to  probe  stability,  a  <  +3%  (of  the  maximum  response  amplitude  in  the  lattice) 
random  kick  was  applied  to  the  lattice  after  each  increment. 

The  upper  boundary  from  o>=.76  to  .89  was  characterized  by  the  lattice  "going  over  the 
top"  after  perturbation  (random  kick).  Because,  the  potential  well  for  the  lower  cutoff  mode 
(Fig.  III.B.3)  is  not  infinitely  high,  the  lattice  is  able  to  escape  the  well  (go  over  the  top)  and 
become  unstable.  There  were  no  visual  indications  in  the  lattice  of  an  impending  instability 
before  the  boundary  was  reached.  The  right  upper  boundary,  from  co  =  .89  to  .94,  was  marked 
by  the  generation  of  another  symmetric  breather  (out  of  phase  with  the  first)  on  the  opposite  side 
of  the  lattice  ring  from  the  initial  breather.  The  lattice  then  went  over  the  top  from  a  random  site 
fairly  quickly  after  the  generated  breather  reached  an  amplitude  comparable  to  the  original 
breather.  Further  tests  with  a  lattice  of  100  sites  indicated  that  the  generated  breather  is 
approximately  27  sites  away  from  the  center  of  the  original  breather  and  will  arise  on  both  sides 
of  the  original  if  the  sites  are  available.  The  breather  structure  in  the  extreme  lower  right  corner 
is  quasiper iodic  until  the  boundaries  where  it  goes  over  the  top  without  the  generation  of  another 
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breather.  The  lower  boundary  from  co=.85  to  .947  is  preceded  by  a  quasiperiodic  region. 
Ultimately  the  lattice  either  goes  over  the  top  or  decays  to  rest  after  very  strong  transients. 
Careful  observations  are  suggestive  of  quasiperiodicity  resulting  from  the  energy  loss  due  to 
soliton  shedding.  The  left  lower  boundary  from  w  =  .76  to  .85  was  only  preceded  by  the 
quasiperiodic  region  at  the  points  indicated  but  the  strong  transients  were  consistently  observed 
before  the  lattice  went  over  the  top.  A  comparison  of  the  stable  region  predicted  by  theory  for 
the  lower  cutoff  symmetric  breather  shows  good  agreement  for  the  upper  right  boundary  but  a 
discrepancy  on  the  lower  boundary  (Fig.  III.B.4).  Some  destabilizing  factor  has  prevented  the 
lower  boundary  from  reaching  the  theoretical  limit.  Experimentation  with  a  half-site  centered 
symmetric  breather  (Fig.  III.B.5)  revealed  a  link  between  the  strong  transients  observed  in  the 
on-site  breather  drive  plane  but  did  not  extend  the  lower  boundary  to  the  theoretical  limit.  The 
strong  transients  occur  at  the  same  drive  amplitudes  that  the  half-site  breather  shows  quasiperiodic 
behavior.  Increasing  the  drive  amplitude  results  in  a  transition  from  half-site  to  on-site  at 
approximately  the  start  of  the  quasiperiodic  line  in  the  on-site  drive  plane.  The  half-site  breather 
is  stable  as  drive  amplitude  is  decreased,  but  it  decays  to  rest  at  approximately  f\ =0.9  for  oj=0.9. 
The  transition  from  steady  state  motion  to  rest  is  interesting  in  that  it  shifts  back  to  an  on-site 
breather  before  it  decays  to  rest. 

Comparison  of  the  theoretical  profile  for  the  parameters  used  in  the  drive  plane  was  done 
for  the  lowest  and  highest  amplitude  saved  profiles  as  well  as  a  half-site  profile.  The  lowest 
amplitude  comparison  was  sharper  and  lower  in  amplitude  than  the  predicted  profile  (Fig. 
III.B.6).  The  highest  amplitude  comparison  was  also  sharper  and  even  lower  in  comparable 
amplitude  (Fig.  II1.B.7).  A  comparison  of  the  half-site  profile  shows  much  closer  agreement  with 
theory  although  the  data  still  exhibits  a  greater  amplitude  (Fig.  III.B.8).  The  amplitude  of  the 
saved  profile  from  the  program  was  not  exactly  the  maximum  amplitude  of  that  particular  state 
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did  not  precisely  coincide  with  the  peak  response  amplitude.  This  was  corrected  by  using  the 
maximum  coordinate  plus  the  coordinates  for  the  preceding  and  following  time  steps  to  find  the 
maximum  amplitude  by  parabolic  curve  fitting.  By  recording  the  time  step,  the  data  could  be 
scaled  according  to  the  time  difference  between  the  turning  point  and  actual  capture  time. 
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Figure  III.B.3   Lower  Cutoff  Potential  Well 
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Symmetric  Half-Site  Lower  Cutoff  Breather 
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Figure  III.B.5   Symmetric  Half -Site  Lower  Cutoff  Breather 
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Lower  Cutoff  Breather  Comparison  To  Theory 
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Figure   III.B.7      Lower  Cutoff   Breather  Theory  Comparison 


Louer  Cutoff  Breather  Comparison  To  Theory 

fi  Q 

8.8 

i              i              i              i 

I 

8.7 

- 

: 

c 

8.6 

- 

u 

"%         8.5 

•m 

1"    8.4 

<c 

| 

8.3 

- 

8.2 

- 

8.1 

J  V 

8             18            28            38            48 

58 

Theory :  1  ine 

Data:  »             Lattice  Site        Bta=.lB8   omega= 

.988 

Figure   III.B.8      Lower   Cutoff   Breather  Theory  Comparison 


29 


C.       UPPER  CUTOFF  BREATHERS 

The  structure  chosen  for  the  upper  cutoff  drive  plane  was  the  symmetric  on-site  breather 
(Fig.  III.C.l).  The  parameters  for  the  drive  plane  were  a=-l  (hardening),  7=.l  (weak 
coupling),  j8=.03  (small)  and  a  lattice  length  of  50  sites.  The  points  in  the  drive  plane  (Fig. 
III.C.2)  were  almost  exclusively  obtained  by  incrementing  the  drive  amplitude  by  .001  and 
recording  the  amplitude  at  which  instabilities  were  observed.  After  each  increment,  a  random 
kick  of  <  +3%  (of  the  maximum  lattice  response  amplitude)  was  applied  to  the  lattice. 

The  upper  left  boundary  marks  the  sharp  transition  from  a  breather  structure  to  a  complex 
upper  cutoff  mode  structure.  The  transition  time  was  fairly  short  and  was  characterized  by  the 
growth  of  one  section  of  the  lattice  to  an  amplitude  comparable  to  that  of  the  original  breather. 
The  growing  section  was  in  a  pure  upper  cutoff  mode  and  spread  to  include  the  rest  of  the  lattice 
after  reaching  full  amplitude.  The  area  inside  the  circles  (Q)  was  where  the  lattice  motion  was 
quasiperiodic.  The  boundary  was  difficult  to  ascertain  for  omega  1.25  and  higher  since  local 
regions  appeared  to  exhibit  quasiperiodicity  but  in  fact  merely  had  long  settling  times.  Regions 
marked  H  denote  areas  of  no  lattice  motion.  Inside  the  region  and  below  the  lower  boundary, 
the  lattice  decayed  to  rest  with  minimal  transients.  Region  QS  was  characterized  by  a 
quasiperiodic  shedding  breather.  The  quasiperiodic  breather  was  usually  observed  near  the 
boundaries  of  QC  and  was  similar  to  Fig.  III.C.l.  The  center  site  of  the  breather  was 
quasiperiodic  and  small  solitons  were  ejected  from  the  center  in  both  directions  simultaneously. 
The  quasiperiodic  breather  spontaneously  transitioned  to  a  complex  cyclic  state  at  the  lower 
boundary  of  QC  (asterisks).  The  lattice  motion  in  the  small  asterisked  area  (QC)  was  of  two 
general  types,  a  complex  cyclic  state  or  an  anti-symmetric  half-site  breather.  The  cyclic  state  that 
evolved  from  the  quasiperiodic  breather  appeared  to  be  chaotic  and  then  settled  into  multiple  on- 
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site  breathers.  The  motion  shifted  back  and  forth  between  the  two  states  with  a  fairly  long  cycle. 
If  the  lattice  was  in  this  particular  mode  when  the  boundary  to  H  (right  side  of  QC)  was  reached, 
the  cycle  was  broken  and  the  lattice  settled  into  a  multiple  breather  state.  The  multiple  breathers 
merged  two  by  two  until  a  single  on-site  breather  was  left.  The  single  breather  then  decayed  to 
rest.  The  lattice  also  spontaneously  shifted  from  the  complex  cycle  to  a  half-site  antisymmetric 
breather  (Fig.  III.C.3).  A  half-site  breather  consistently  evolved  from  an  quasiperiodic  on-site 
breather  at  the  upper  boundary  of  QC.  The  half-site  breather  exhibited  a  degeneration  similar 
to  that  of  the  complex  cyclic  at  the  boundary  to  region  H.  The  breather  spontaneously  changed 
back  into  an  on-site  breather  and  then  decayed  to  rest.  Multiple  antisymmetric  half-site  breathers 
were  also  observed.  A  state  with  two  half-sites  10-15  sites  apart  (center-to-center)  were  believed 
to  be  transferring  energy  by  soliton  shedding.  A  growth  in  amplitude  in  the  center  sites  of  one 
breather  corresponded  to  the  arrival  of  a  soliton  from  the  other  breather  and  the  expected 
decrease  in  the  center-site  amplitude  of  the  second  breather  due  to  soliton  ejection  was  observed. 

The  upper  and  lower  boundaries  compare  well  with  the  theoretical  boundaries  (Fig. 
III.C.4).  The  upper  boundary  does  not  quite  follow  the  hyperbola  at  the  higher  drive  amplitudes, 
however  the  theory  is  approximate  at  these  amplitudes.  The  lower  boundary,  which  corresponds 
to  the  lattice  decaying  to  rest,  coincides  almost  exactly  with  the  predicted  threshold  of  lattice 
motion.  Comparison  of  theoretical  modulation  envelopes  and  the  numerical  data  were 
completed  for  extremely  high,  high  and  low  response  amplitude  cases  as  well  as  a  half-site  case. 
The  observed  amplitudes  were  adjusted  to  that  of  the  turning  point  as  described  in  Section  II. B. 
The  normal  180°  phase  difference  between  lattice  sites  was  eliminated  for  easier  modulation 
envelope  comparison.  The  extreme  case  was  not  as  sharp  as  predicted  and  was  =70%  of  the 
theoretical  amplitude  (Fig.  III.C.5).    Similarly,  the  amplitude  of  the  high  case  was  =90%  of 
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theoretical  but  the  envelope  was  sharper  than  expected  (Fig.  III.C.6).  The  low  amplitude  case 
exhibited  the  inverse  characteristics  of  the  extreme  case  (Fig.  III.C.7).  The  observed  envelope 
was  sharper  than  the  theory  and  =  5%  higher  in  amplitude.  A  comparison  of  the  half-site  case 
showed  very  good  agreement  with  theory  (Fig.  III.C.8). 
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D.       SOLITON  SHEDDING 

The  initial  states  for  both  upper  and  lower  cutoff  mode  response  investigations  were 
evanescent  states  in  a  lattice  of  100  sites.  One  end  of  the  lattice  was  driven  and  the  other  end 
was  given  a  free  boundary  condition.  The  large  number  of  sites  minimized  reflection  while 
maintaining  the  speed  of  the  simulation.  The  parameters  for  the  upper  cutoff  were  7=. 75, 
0=.O3,  while  the  lower  cutoff  were  7=. 50,  jS=.03.  Snapshots  in  time  of  the  on-screen  motion 
of  the  lattice  were  obtained  to  help  describe  the  motion.  In  these  "snapshot"  figures,  each  dot 
is  an  individual  lattice  site  and  the  site  furthest  to  the  left  is  the  driven  end  site. 

Soliton  shedding  was  observed  in  both  the  upper  and  lower  cutoff,  end  driven  lattices.  The 
shedding  was  similar  for  both  modes,  with  differences  in  the  actual  method  of  energy  transfer 
from  the  driven  site.  The  upper  cutoff  shedding  was  produced  by  the  cyclic  relaxation  in 
amplitude  of  the  driven  end  site.  The  lower  cutoff  shedding  was  also  generated  by  end  site 
relaxation  but  transferred  energy  through  a  "pivot"  point  before  the  soliton  was  created. 

An  artificial  potential  well  of  the  form: 

fix)    -  ~X  III.D.l 

y/l  +  x2 

was  used  to  replace  f(x)=  -x  +  x3  in  the  lower  cutoff  program  to  permit  amplitudes  of  the 

magnitude  necessary  to  observe  shedding  in  the  softening  mode.  Otherwise,  the  oscillators  would 

"go  over  the  top."    The  upper  threshold  was  marked  by  a  distinct  motion  change  and  soliton 

shedding  (Fig.  III.D.l).    The  shed  solitons  were  very  small  and  strongly  evanescent  from  a>=  .95 

down  to  a>=  .6.   The  soliton  shedding  at  o>=  .6  was  distinct  and  the  structure  easily  visible  until 

20-30  sites  from  the  driven  end.  The  shedding  results  were  similar  for  co  =  .5  with  solitons  visible 

until  approximately  site  35.    A  third  harmonic  component  appeared  at  low  drive  amplitudes 
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(r/= .  144)  and  continued  until  the  shedding  threshold.  The  harmonic  was  identified  by  taking  the 
FFT  of  a  time  series  from  the  driven  site  (Fig.  III.D.2).  The  visible  result  of  the  harmonic  was 
an  evanescent  ripple  wave  (Fig.  III.D.3). 

A  fifth  harmonic  was  identified  for  co=.3  (Fig.  III.D.4)  and  it  also  disappeared  at  the 
shedding  threshold.  The  strength  of  the  harmonic  appeared  to  vary  slowly  with  amplitude  but 
over  a  fairly  large  range.  The  driven  end  site  became  quasiperiodic  from  77  =  .415  to  1.42  and 
from  2.14  to  the  shedding  threshold  for  co  =  .2.  Small  wave  packet  ejection  was  observed  from 
tj=2.00  to  2.59.  Upon  reaching  the  shedding  threshold,  the  end  site  relaxed  and  a  large  soliton 
was  shed  (Fig.  III.D.5).  The  soliton  would  propagate  6-8  sites,  then  appear  to  shed  a  much 
smaller  soliton.  Both  the  original  and  shed  soliton  would  then  damp  out  quickly  (Fig.  III.D.6). 
A  seventh  harmonic  appeared  for  to  =  .2  but  was  not  visually  identifiable.  Odd  harmonics  up  to 
the  11th  were  identified  for  w=.l  (Fig.  III. D. 7-8)  and  visually  affected  the  lattice  motion.  The 
end  site  quasiperiodicity  started  at  77  =  2.25  and  the  evanescent  wave  packets  started  at  77=  1.35. 
The  shedding  threshold  was  not  reached  by  17=6.0. 

The  dynamic  threshold  for  the  upper  cutoff  lattice  exhibits  a  consistent  trend  even  though 
the  motion  of  the  resultant  state  changes  significantly  between  tj  =  2.32  and  2.33  (Fig.  III.D.9). 
The  basic  motion  changes  from  soliton  shedding  to  a  non-shedding  evanescent  state.  Shedding 
for  co =2.025  occurs  very  slowly  and  the  soliton  propagates  30-40  lattice  sites  and  then  appears 
to  stop  (Fig.  III.D.  10).  Apparently  the  nonlinearity  lowers  the  frequency  enough  that  linear  wave 
propagation  can  occur  inside  the  structure  and  it  appears  to  "ring"  as  it  decays  (Fig.  III.D.  11). 
The  shedding  occurs  slowly  for  o>=2.05  and  the  soliton  is  large  but  jumbled  (Fig.  III.D.  12)  and 
dies  out  quickly.  The  size  of  the  shed  solitons  is  cyclic  for  co=2. 1-2.3  and  the  structure  is  very 
diffuse  (Fig.  III.D.  13).  Several  small  solitons  are  ejected  followed  by  one  large  one  and  the 
relaxation  of  the  end  site  is  very  distinct.  Intermittent  shedding  occurs  for  w=2.31  and  the  first 
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8-9  sites  move  together  (Fig.  III. D.  14).  The  transition  point  from  shedding  to  nonshedding 
seems  to  be  between  co=2.32  and  2.33.  The  motion  for  o>  =  2.32  is  varied  in  the  number  of  sites 
in  tandem,  the  size  of  the  shed  solitons  and  the  magnitude  of  the  end  site  relaxation.  The 
shedding  is  transitional  for  w=2.33  and  the  motion  settles  into  a  complex  evanescent  state.  The 
first  6-9  sites  move  in  tandem  and  the  last  site  is  the  first  of  2-3  evanescent  sites.  The  motion 
for  w  =  2.35-2.6  is  similar  but  only  2-3  sites  are  in  tandem.  The  vertical  boundary  between 
shedding  and  non-shedding  states  was  confirmed  by  obtaining  high  drive  level  evanescent  states 
and  crossing  the  boundary  into  the  shedding  region.  Specifically,  the  states  w=2.35  77=5.05  and 
w  =  2.5  t/  =  8.0  still  exhibited  the  same  lattice  motion  initiatated  at  the  transition  boundary.  The 
states  were  then  decremented  in  frequency  to  cross  into  the  shedding  region.  The  motion  of  both 
states  slowly  evolved  to  the  previously  described  lattice  motion  for  the  final  frequency. 
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Figure   III.D.4      Lower  Cutoff   FFT   Spectrum  w=.3 
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Figure   III.D.5      Time   Snapshots   co=.2 
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Figure  III.D.6   Time  Snapshots  w=.2 
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Figure   III. D. 10      Time   Snapshots   w=2.025 
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Figure    III. D.  11      Time   Snapshot   co=2.025 
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Figure  III. D.  12   Time  Snapshot  to=2.05 


Eta-2  2   Omeg«-2  2 


Figure  III. D. 13   Time  Snapshot  w=2 . 2 
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Figure  III. D.  14   Time  Snapshot  o>=2.31 


Figure  III. D.  15   Time  Snapshot  co=2.4 
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IV.   SUMMARY  AND  CONCLUSIONS 

It  has  been  shown  numerically  that  nonpropagating  steady  state  breather  solitons  can  occur 
in  one-dimensional  lattices  of  coupled  nonlinear  oscillators  that  are  damped  and  parametrically 
driven.  The  simplicity  of  the  model  suggests  that  these  states  can  exist  in  a  variety  of  systems. 
Amplitude  data  agree  well  with  a  nonlinear  Schrodinger  theory  at  low  amplitudes,  but  disagree 
substantially  at  higher  amplitudes.  The  theory  assumes  that  the  amplitudes  are  weakly  nonlinear 
and  slowly  varying  in  space.  The  continued  existence  of  breathers  at  amplitudes  where  the  theory 
breaks  down  is  evidence  of  the  robustness  of  these  states.  Dissipation  and  drive  stabilize  the 
breathers,  although  this  is  not  yet  understood. 

Breathers  exist  for  a  range  of  drive  parameters  (amplitude  and  frequency),  and  these 
regions  have  been  mapped  for  both  lower  and  upper  cutoff  breathers.  Each  region  substantially 
disagreed  with  the  theoretical  prediction  and,  moreover,  the  natures  of  the  two  numerical  regions 
were  strongly  dissimilar.  The  lower  cutoff  breathers  are  found  to  have  a  low-amplitude 
instability  that  is  not  predicted  by  the  theory.  The  drive  parameter  region  of  the  upper  cutoff 
breathers  contains  a  relatively  large  "finger"  in  which  quasiperiodicity  occurs,  and  a  relatively 
large  "island"  in  which  the  motion  decays  to  rest.  The  corresponding  instabilities  are  not  yet 
understood,  but  it  appears  that  the  quasiperiodicity  is  a  consequence  of  soliton  shedding  even 
though  the  shedding  is  not  apparent  except  in  a  small  region  of  the  drive  plane.  This  is  based 
on  the  observation  of  strong  quasiperiodic  amplitude  modulation  during  shedding  and  relatively 
weak  modulation  when  shedding  is  not  visible.  Clearly  visible  soliton  shedding  is  not  expected 
to  occur  for  a  global  parametric  drive  because  a  propagating  breather  has  a  spatially  varying 
phase,  whereas  the  drive  is  spatially  constant. 
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For  an  end  driven  lattice,  soliton  shedding  is  dramatically  visible.  There  is  a  drive 
amplitude  threshold  for  the  shedding  to  occur,  and  the  value  decreases  as  the  linear  frequency  is 
approached.  The  phenomenon,  although  not  yet  understood,  is  expected  to  occur  in  any 
nonlinear  wave  system  that  possesses  breathers.  A  possible  application  is  the  regeneration  of 
fiber  optic  solitons  which  attenuate  with  distance. 
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APPENDIX 

A.       PROGRAM  COMMANDS 

1.  a  -  Manual  frequency  increment,  adjustable  with  A. 

2.  A  -  Manual  frequency  increment  adjustment,  can  be  positive  or  negative. 

3.  d  -  Coarse  damping  adjustment  (decreasing),  nominally  .01,  can  be  changed  in 
variable  declaration  section  of  source  code  (the  variable  is  Betaincrement). 

4.  D  -  Coarse  frequency  adjustment  (decreasing),  nominally  .001,  can  be  changed  in 
variable  declaration  section  of  source  code  (Omegaincrement). 

5.  Ctrl  D  -  Coarse  drive  amplitude  adjustment  (decreasing),  nominally  .001,  can  be 
changed  in  variable  declaration  section  of  source  code  (Etaincrement). 

6.  e  -  Fine  damping  adjustment  (decreasing),  1/10  of  coarse  damping  increment,  can 
be  changed  in  the  mainO  section  of  the  source  code  under  the  appropriate  key  hit 
section. 

7.  E  -  Fine  frequency  adjustment  (decreasing),  1/10  of  coarse  frequency  increment,  can 
be  changed  in  the  mainO  part  of  the  source  code  under  the  appropriate  key  hit 
section. 

8.  Ctrl  E  -  Fine  drive  amplitude  adjustment  (decreasing),  1/10  of  coarse  drive 
amplitude  increment,  can  be  changed  in  the  mainO  section  of  the  source  code  under 
the  appropriate  key  hit  section. 

9.  Ctrl  F  -  Time  series  record  of  one  lattice  site,  prompts  for  a  file  name  to  store  the 
series.     This  function  works  only  in  Text  Mode  (Ctrl  T).     The  number  and 
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periodicity  of  the  samples  are  controlled  in  the  eqnmotionO  function  section  of  the 
source  code.   Set  for  80  samples,  taken  at  the  rate  of  one  per  time  increment. 

10.  G  -  Coupling  adjustment  (gamma),  displays  current  value  asks  for  new  value. 

1 1 .  Ctrl  G  -  Graphics  mode,  displays  real  time  lattice  motion.  Also  refreshes  the  screen 
without  interfering  with  the  lattice  motion.  Automatically  invoked  when  the 
program  is  started  or  restarted. 

12.  Ctrl  H  -  Total  lattice  energy  monitor,  can  only  be  used  in  Text  Mode. 

13.  i  -  Zoom  in.  The  amplification  of  the  screen  is  2n  for  n  button  pushes.  The 
amplification  is  vertical  only. 

14.  Ctrl  I  -  Increases  the  strobing  frequency  of  the  phaseplot  option,  effectively 
increasing  the  sampling  rate. 

15.  Alt  K  -  Coupling  modulation,   Options  are  random,  gradient  or  sine  wave. 

16.  Ctrl  K  -  Amplitude  kick,  specified  lattice  site  by  the  desired  amount.  Amplitude  and 
lattice  site  input  from  the  keyboard. 

17.  Ctrl  L  -  Pins  specified  lattice  site.  Acts  like  a  toggle  switch,  reselection  unpins  the 
site. 

18.  n  -  Random  perturbation  of  the  lattice.  Perturbs  all  sites  with  a  random  coordinate 
change  of  <  +  3%  of  the  maximum  amplitude  in  the  lattice. 

19.  o  -  Zoom  out.  The  reduction  of  the  screen  is  2n  for  n  button  pushes.  The  reduction 
is  vertical  only. 

20.  Ctrl  O  -  Decreases  the  strobe  frequency  of  the  phaseplot  option,  effectively 
decreasing  the  sampling  rate. 
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21.  p  -  Pause  program.  Freezes  execution  of  the  program.  Alternate  method  of 
obtaining  a  time  series  file. 

22.  P  -  Phaseplot  mode.    Selects  up  to  five  lattice  site  for  phase  space  monitoring. 

23.  Ctrl  Q  -  Exits  the  program. 

24.  r  -  User  initiated  state  save.  The  lattice  must  be  stable  to  within  1  %  (blue  color)  or 
the  option  will  freeze  the  program  until  the  1  %  criterion  is  met.  The  option  saves 
the  current  lattice  modulation  at  the  upper  turning  point  of  the  lattice  site  being 
monitored. 

25.  Ctrl  R  -  Restarts  the  program  without  going  back  out  to  DOS.  Goes  to  the  same 
screen  that  was  displayed  when  the  program  was  first  started  up. 

26.  s  -  Program  freeze,  captures  lattice  modulation  when  pushed.  Any  key  will  continue 
the  program  but  the  lattice  will  be  at  rest. 

27.  S  -  Monitors  the  stability  of  the  specified  lattice  site.  Works  like  a  toggle  switch. 
The  lattice  will  change  colors  when  the  monitored  site  is  within  1  %  of  its  average 
amplitude  for  the  last  twenty  upper  turning  points. 

28.  Ctrl  S  -  Turns  off  the  drive. 

29.  Ctrl  T  -  Text  Mode.  Prints  system  parameters  on  screen.  Total  lattice  energy  can 
be  monitored  in  this  mode. 

30.  u  -  Coarse  damping  adjustment  (increasing),  nominally  .01,  can  be  changed  in 
variable  declaration  section  of  source  code  (variable  is  Betaincrement). 

31.  Ctrl  U  -  Coarse  drive  amplitude  adjustment  (increasing),  nominally  .001,  can  be 
changed  in  variable  declaration  section  of  source  code  (Etaincrement). 
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32.  U  -  Coarse  Frequency  adjustment  (increasing),  nominally  .001,  can  be  changed  in 
variable  declaration  section  of  source  code  (Omegaincrement). 

33.  v  -  Fine  damping  adjustment  (increasing),  1/10  of  coarse  damping  increment,  can 
be  changed  in  the  mainO  section  of  the  source  code  under  the  appropriate  key  hit 
section. 

34.  Ctrl  V  -  Fine  drive  amplitude  adjustment  (increasing),  1/10  of  coarse  frequency 
increment,  can  be  changed  in  the  mainO  section  of  the  source  code  under  the 
appropriate  key  hit  section. 

35.  V  -  Fine  drive  amplitude  adjustment  (increasing),  1/10  of  coarse  increment,  can  be 
changed  in  the  mainO  part  of  the  s     rce  code  under  the  appropriate  key  hit  section. 

36.  w  -  Displays  peak  lattice  amplitude,  works  in  Graphics  Mode  only. 

37.  Ctrl  W  -  Turns  on  waterfall  type  display.  Trends  are  easily  noticed  on  this  type  of 
display. 

38.  Ctrl  X  -  Selects  lattice  site  to  monitor  for  spectrum  analysis. 

39.  y  -  Dump  spectrum  generated  by  Ctrl  X. 

40.  z  -  Sets  damping  to  zero. 

41.  Z  -  Sets  frequency  to  zero. 

42.  Ctrl  Z  -  Zeros  the  peak  amplitude  variable  (variable  automatically  zeros  when  a 
parameter  is  changed). 
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B.       PROGRAM  CODE 

The  QuickC  program  codes  for  the  aforementioned  numerical  implementations  (Sect  III) 
are  listed  below.  The  full  program  listing  for  Lower  Cutoff  (CCLATL.C)  is  included  and  the 
equation  of  motion  for  Upper  Cutoff  (CCLATU.C),  Lower  Cutoff  end-driven  (ENDL.C)  and 
Upper  Cutoff  end-driven  (ENDU.C). 

1.        Lower  Cutoff,  Global  Drive 

/*      PROGRAM  LATTICE  (VGA) 
VERSION  2.0  (QUICKC) 
*/ 

#define  LASTUPDATE  22  JUL  1991  BY  CLEON  WALDEN 
/* 

THIS  PROGRAM  SIMULATES  A  GENERAL  LATTICE  WITH  EQUATIONS  OF  MOTION 
THAT  CAN  BE  SUBSTITUTED  IN  WHERE  INDICATED.    VARIOUS  INTERACTIVE 
FEATURES  ARE  PROVIDED,  WHICH  ARE  EXPLAINED  IN  THE  PROGRAMMER'S 
MANUAL  FOR  LATTICE.  ENERGY  MONITORING  IS  ADDED  TO  TEXT  SCREEN. 
A  WATERFALL  DISPLAY  IS  ADDED  TO  MONITOR  SOLITON  MOTION  DUE  TO 
MEDIUM  EFFECTS. 
*/ 

#include  "BIOS.H" 
#include  "graph. h" 
#include  "stdlib.h" 
include  "CONIO.H" 
include  "MATH.H" 
include  "STDIO.H" 
#include  "time.h" 
^include  "direct.h" 
^include  "process. h" 
^include  "dos.h" 

#define  getrandom(min,max)  ((randO  %  (int)((max)-(min)))+(min)+l) 

#define  SQR(a)  ((a)*(a)) 

#define  CUB(a)  ((a)*(a)*(a)) 

#define  SWAP(a,b)  tempr=(a);(a)  =  (b);(b)=tempr 

#define  DOFOR(i,to)  for(i=0;i<to;i+  +) 

#define  PI  3.14159265359 

#define  SCREENCORRECTIONFACTOR  1.05 

#define  etaincrement  0.01 

#define  beta  increment  0.01 
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#define  omegaincrement  0.001 

#define  STABILITYINCREMENT  0.01 

#define  DISPLAYINCREMENT  DINC 

#define  PHASEJNCREMENT  PINC 

#define  MIDDLE  ELEMENT  (no_pendulums/2) 

#define  textflag  flags[0] 

#define  graphicsflag  flags[l] 

#define  filedumpflag  flags[2] 

#define  stability  flag  flags[3] 

#define  peakflag  flags  [4] 

#define  stableflag  flags[5] 

#define  phaseflag  flags[6] 

#define  pauseflag  flags  [7] 

#define  pause_flag2  flags[8] 

#define  stopflag  flags[9] 

#define  spectrumflag  flags[10] 

#define  dumpspectrumflag  flagsfll] 

#define  energyflag  flags[12] 

#define  waterfallflag  flags[13] 

^define  wait_flag  flags[14] 

/*      The  dynamical  variables  are  entered  here.    */ 

double  coordinate!  1 50] , momentum [  1 50] ,old_coordinate[  1 50] ,old_momentum[  150] , 
oldold_coordinate[  1 50]  ,oldold_momentum[  1 50] ; 

double  acceleration,gamma[  150],mean_gamma,beta,omega0[  150],mean_omega0, 
eta,omega,alpha,model_time,time_int, 

max_amp,mode_amp,time_series_disp[8000],phase_elements[5], 
peak_record[20],pinned_elements[150],DINC,PINC,period, 
spectrum[4098],temp2[4096],energy,gradient,wat_inc, 
peak_amp,amp_kick,freq_inc,delta; 

int     no_pendulums,flags[  15], counter  l,phase_counter,phase_int,first_element, 
chosen_element,stability_element,disp_color,colors[5],counter2,step_size, 
spectrum_element,spectrum_counter,sample_counter,energy_counter,g,gst[8000], 
waterfall_counter,color_counter,number_click:s,p_flag,t_flag,top_flag; 

mainO  { 

FILE  *fp; 

int  c,ij,k,l,m,n,xx,a; 

double  modulation; 

char  ans[5],buffl50],buff2[50]; 

/*    void  display _textO,stability_restartO,user_initO,display_graphicsO, 

process_pauseO,monitor_stabilityO,start_file_dumpO, 

Stop_file_dump0,pin_elements0,kick_elements0,stability_check0, 

init_phasplotO,phasplotO,compute_spectrumO,plot_spectrumO;  */ 

_clearscreen(wGCLEARSCREEN); 


56 


userinitO; 
_setvideomodeCMRES4COLOR); 

setcolor(l); 
sprintf(buff,"Eta  %.31f  Omega  %.31f  Gamma  %.21f',eta,omega,mean_gamma); 

_registerfonts("TMSRB.FON"); 

_moveto(10,10); 

_setfont("t'tms  rmn'bn5"); 

_outgtext(buff); 
xx =0; 
wat_inc=l; 
disp_color=l; 
stop_flag=0; 
counter2=0; 
energy_flag=0; 
energycounter = 0; 
color_counter=0; 
while(stop_flag==0)  { 
energy =0; 

if(pause_flag2=  =  1)  { 
pause_flag2=0; 
pause_flag=0; 
process_pauseO; 
}   /*  if(pause...)  */ 
if(kbhit0!=0){ 
c=getch0; 
if(c=  =  17){ 
stop_flag=l; 

} 

else  if(c=  =  112)  { 

pause_flag=l; 

} 

else  if(c==20)   { 

_clearscreen(_GCLEARSCREEN); 

text_flag=l; 

graphics  flag =0; 

waterfall_flag=0; 

wait_flag=0; 

spectrum_flag=0; 

displaytextO; 

} 

elseif(c==7)   { 

screengraphO; 

phase  flag =0; 

wait_flag=0; 

spectrum_flag=0; 

waterfall_flag=0; 

text_flag=0; 


/: 


/*  Check  for  keystroke  from  user  */ 
AQ  :   Quit  program  */ 


/*p 


/" 


Pause  program/dump  state  */ 


Text  Mode  */ 


*  *t 


Graphics  Mode  */ 
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graphicsflag  =  1 ; 

sprintf(buff,"Eta  %.31f  Omega  %.31f  Gamma  %.21f",eta,omega,mean_gamma); 

_registerfonts("TMSRB.FON"); 

_moveto(10,10); 

_setfont("t'tms  rmn'bn5"); 

_outgtext(buff); 

} 

elseif(c==8)   {       /*   AH  :    Monitor  energy  in  text  mode  */ 

if(energy_flag=  =0)  energyflag  =  1 ; 

el  se  energyfl  ag = 0 ; 

} 

else  if (c  =  =  16)   {      /*  AP         :   Kick  in  parameter  variations  */ 

screentextO; 

srand((unsigned)time(NULL)); 
printf("\nEnter  1  if  you  wish  to  vary  coupling:    "); 
scanf("%d",&j); 
if(j=  =  l){ 
printf("\nEnter  1  (random),  2  (gradient),  or  3(sine)  desired:    "); 
scanf("%d",&k); 
if(k==2)  { 
printf("\nEnter  gradient  (percentage  over  entire  lattice):    "); 
scanf("%lf",&gradient); 
DOFOR(i,no_pendulums)  { 

gamma[i] =gamma[i]  +  mean_gamma*i*gradient/(100*no_pendulums); 
}  /*  DOFOR  */ 
}  /*  if(j  =  =2)  */ 
elseif(k==3)  { 

printf("\nEnter  modulation  amplitude  (percentage):  "); 
scanf("%lf',&gradient); 

printf("\nEnter  integer  number  of  wavelengths  to  use:    "); 
scanf("%d",&l); 
DOFOR(i,no_pendulums)   { 

gamma[i]=gamma[i]  +  mean_gamma*gradient/100*(-cos(2*l*PI*i/ 
no_pendulums)); 

}" 
} 
else  DOFOR(i,no_pendulums)  { 

j=randO; 

gamma[i] =mean_gamma*(.95  + .  1  *j/RAND_MAX); 

} 
} 
else  { 

printf("\nVarying  omegao...\n"); 

printf("\nEnter  1  (random),  2  (gradient),  or  3(sine)  desired:    "); 

scanf("%d",&k); 

if(k==2){ 
printf("\nEnter  gradient  (percentage  over  entire  lattice):    "); 


58 


scanf("%lf',&gradient); 
DOFOR(i,no_pendulums)  { 
omegaO[i]=omegaO[i]  +  mean_omegaO*i*gradient/(100*no_pendulums); 

} 

} 

elseif(k==3)  { 
printf("\nEnter  modulation  amplitude  (percentage):    "); 
scanf("%lf*,&gradient); 

printf("\nEnter  integer  number  of  wavelengths  to  use:    "); 
scanf("%d",&l); 
DOFOR(i,no_pendulums)  { 
omegaO[i]=omegaO[i]+mean_omegaO*gradient/100*(-cos(2*l*PI*i/no_pendulums)); 

} 

} 
else  { 

DOFOR(i,no_pendulums)  { 
j  =  randO; 

printff  Random  number  is:  %d  RAND_MAX  is  %d\n",j,RAND_MAX); 
omegaO[i]  =  mean_omega0%95  + .  1  *j/RAND_M  AX); 

} 
} 
} 
screen_graphO; 

} 

elseif(c==23)   {       /*  "W  : 

if(waterfall_flag==0)  { 

phase_flag=0; 

graphics_flag=0; 

text_flag=0; 

spectrum  flag =0; 

waterfall_flag=l; 

initwaterfallO; 

} 
else  { 

if(wait_flag=  =  l)  { 

wait_flag=0; 

initwaterfallO; 

} 

else  waterfall_flag=0; 

} 
} 
elseif(c==80)   {        /*  P  : 

if(phase_flag==0)  init_phasplotO; 

else  { 

phase_flag=0; 

spectrum_flag=0; 

graphics_flag=l; 


Waterfall  Display  Mode  */ 


Phase  Plot  Mode  */ 
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waterfall_flag=0; 
wait_flag=0; 
screengraphO; 
}  /*  else  */ 

} 

elseif(c==6)   {  /*  AF         :   Time  series  to  file  */ 

if(file_dump_flag=  =0)  startfiledumpO; 
else  { 

stopfiledumpO; 
filedumpflag = 0; 

} 
} 

else  if(c==47)   {  /*  G         :   Change  gamma  */ 

screentextO; 

printf("\nEnter  new  gamma,  gamma  is  %lf',mean_gamma); 
scanf("  %  If '  ,&mean_gamma); 
DOFOR(i,no_pendulums)  { 
gamma[i]  =  meangamma; 

} 

screengraphO; 

stabilityrestartO; 

} 

else  if (c  =  =  119)   {  /*  w         :   Peak  amplitude  */ 

_moveto(210,10); 
_setfont("t'tms  rmn'bn5"); 
sprintf (buffi, "Amp  %  .41f ' ,peak_amp); 
_outgtext(buff2); 

} 

else  if(c==26)   {  /*  AZ         :   Reset  peak_amp  variable  */ 

peak_amp=0.0; 
} 


else  if(c  =  =  24)   {        /*  AX         :   Calculate  spectrum  */ 

spectrum  _flag=l; 

screen_textO; 

printf("\nEnter  number  of  element  to  be  analyzed:    "); 

scanf("  %  d "  ,&spectrum_element); 

screengraphO; 

} 

else  if (c  ==21)   {        /*  AU  :   Increase  drive  (coarse)  */ 

savestateO; 
eta = eta + etaincrement; 
stabilityrestartO; 

} 

elseif(c==4)   {  /*  AD         :   Decrease  drive  (coarse)*/ 

save_stateO; 
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eta = eta-etaincrement; 
stabilityrestartO; 

} 

else  if(c==22)   {        /*  AV  :   Increase  drive  (fine)   */ 

savestateO; 

eta = eta + etaincrement/ 10; 
stabilityrestartO; 

} 

else  if(c==5)   {         /*  AE         :   Decrease  drive    (fine)*/ 

save_stateO; 

eta = eta-etaincrement/ 10; 
stabilityrestartO; 

} 

elseif(c==26)    {        /*  *S  :   Turn  off  drive  */ 

savestateO; 
eta=0; 

stabilityrestartO; 

} 

else  if(c=  =85)   {        /*  U  :    Increase  frequency  (coarse)*/ 

save_stateO; 
timewaitO; 

time_int=time_int*200/period; 
omega = omega  -I-  omegaincrement; 
period  =  2*PI/omega; 
timeint = time_int*period/200; 
stabilityrestartO; 

} 

else  if(c==68)   {        /*  D  :  Decrease  frequency  (coarse)*/ 

savestateO; 
timewaitO; 

timeint = time_int*200/period; 
omega = omega-omegaincrement; 
period = 2  *PI/omega; 
timeint = time_int*period/200; 
stability_restartO; 

} 

else  if(c  =  =  86)   {        /*  V  :   Increase  frequency  (fine)*/ 

save_stateO; 
timewaitO; 

time_int=time_int*200/period; 
omega = omega  -I-  omegaincrement/ 1 0; 
period = 2  *PI/omega; 
timeint = time_int*period/200; 
stabilityrestartO; 

} 

else  if(c==69)   {        /*  E  :   Decrease  frequency  (fine)*/ 

savestateO; 
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/*z 


Set  frequency  to  zero  */ 


:   Set  manual  frequency  increment*/ 


Manual  Frequency  increment*/ 


timewaitO; 

timeint = time_int*200/period ; 

omega = omega-omegaincrement/ 1 0; 

period =2  *PI/omega; 

timeint = time_int*period/200; 

stabilityrestartO; 

} 

elseif(c==90)   { 

savestateO; 
omega =0; 
stabilityrestartO; 

} 

elseif(c==65)     {       /*  AA 

screentextO; 

printf(" Manual  frequency  increment  is  %lf\n",freq_inc); 

printf("New  increment:"); 

scanf("%lf',&freq_inc); 

screen_graphO; 

} 

else  if(c==97)     {       /*  a 

save_state(); 
timewaitO; 

timeint = timeint  *200/period ; 
omega = omega + freq_inc; 
period  =  2*Pl/omega; 
timeint = time_int*period/200; 
stabilityrestartO; 

} 

else  if(c=  =  117)   {       /*  u 

save_stateO; 
beta = beta + betaincrement ; 
stabilityrestartO; 

} 

else  if(c=  =  100)   {        /*  d 

savestateO; 
beta = beta-betaincrement; 
stability_restartO; 

} 

else  if(c=  =  118)   {         /*  v 

savestateO; 
beta = beta + betaincrement/ 10; 
stabilityrestartO; 

} 

else  if(c=  =  121)  {         /*  y 

if  (spectrum  flag  =  =  1)  dump_spectrum_flag=  1; 

} 

else  if(c=  =  101)   {         /*  e 


Increase  damping  (coarse)*/ 


Decrease  damping  (coarse)*/ 


Increase  damping  (fine)*/ 


Dump  spectrum   */ 


:   Decrease  damping  (fine)*/ 
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Turn  off  damping  */ 


Pin  elements  */ 


Kick  elements  */ 


savestateO; 
beta = beta-betaincrement/ 1 0; 
stabilityrestartO; 

} 

else  if(c=  =  122)   {  /*  z 

save_state(); 
beta=0; 
stabilityrestartO; 

} 

elseif(c=  =  12)   {  /*  "L 

pinelementsO; 

stabilityrestartO; 

} 

else  if(c=  =  ll)   {  /*  AK 

kickelementsO; 
stabilityrestartO; 

} 

else  if(c=  =  105)  {  /*  i  :   Zoom  in   */ 

DINC  =  DINC/2; 

PINC  =  PINC*2; 

watinc = watinc  *2 ; 

} 

else  if(c=  =  1 1 1)  {  /*  o  :   Zoom  out  */ 

DINC=DINC*2; 
PINC  =  PINC/2; 
watinc  =  wat_inc/2; 

} 

else  if(c=  =  110)  {  /*  n 

/*  Maximum  +/-3%  of  peak  value  */ 

linearkickO; 

screen_graphO; 

stabilityrestartO; 

} 

elseif(c=  =  15)   {  /*  A0 

phaseint = phaseint- 1 ; 

} 

elseif(c==9)     {  /*  AI 

phaseint = phaseint  + 1 ; 

} 

else  if(c=  =  18)   {  /*  "R 

screentextO; 

userinitO; 

sprintf(buff,"Eta  %.31f  Omega  %.31f  Gamma  %.21f',eta,omega,mean_gamma); 

_registerfonts("TMSRB.FON"); 

_moveto(10,10); 

_setfont("t'tms  rmn'bn5"); 

_outgtext(buff); 


Perturb  system  */ 


:   Decrease  strobe  freq  */ 
Increase  strobe  freq  */ 
Restart  program  */ 
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if(phase_flag=  =0)  _setvideomode(_MRES4COLOR); 
else  _setvideomodeCVRES16COLOR); 

} 

else  if(c=  =  1 15)   {  /*  s  :   Stop  everything  as  is  */ 

DOFOR(i,no_pendulums)  { 

coordinate[i]=0; 

momentum[i]=0; 

}  /*  DOFOR  */ 

stabilityrestartO; 

} 

else  if(c==83)  {  /*  S  :   Monitor  stability  */ 

monitorstabilityO; 

sprintf(buff,"Eta  %.31f  Omega  %.31f  Gamma  %.21f",eta,omega,mean_gamma); 

_registerfonts("TMSRB.FON"); 

_moveto(10,10); 

_setfont("t'tms  rmn'bn5"); 

_outgtext(buff); 

} 

else  if(c=  =  114)   {   /*   r         :    User  initiated  state  save  */ 

save_stateu(); 

} 

}      /*  if(kbhit...)  */ 

if(wait_flag==0)  { 

eqnmotionO; 

}  /*  if(wait  flag  */ 

}  /*  while(stop_flag...  */ 

_setvideomode(_DEFAULTMODE); 

printf("PROGRAM  COMPLETE  AT  %lf,model_time); 

}    /*  MAINO  */ 

userinitO  { 

FILE  *fq; 

int  c,i,j,k,l,m; 

char  answer[  1  ]  ,ans[  1  ]  ,answ[  1  ], filename! 30] ; 

printf(" GENERALIZED  LATTICE  MODEL  PROGRAM  W/  VARIABLE  PARAMETERS^"); 

printf(" Version  2.1X486  (MS)  \n"); 

printf("Last  updated  22  JUL  1991\n"); 

printf("Variant  notes:  Default  IC  is  AM\n"); 

printf("Omega0  is  set  equal  to  one  for  all  cases!\n"); 

printf(" Rotating  phase  plane  is  usedAn"); 

printf("Real  time  FFT  function  is  added... \n"); 

printf(" Energy  monitoring  available  via  AH  in  text  mode\n"); 

printf("Set  topflag  to  get  output  files  for  theory  comparison. \n"); 

printf(" Manual  Frequency  increment  is  set  to  1/20  omega  inc.\n"); 

/*  Sets  flag  to  output  additional  information  for  theory  comparison. 

Adds  the  information  to  the  end  of  the  normal  output  file*/ 
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printf("Set  top_flag?   "); 

scanf("%s",ans); 

if((ans[0]  =  =  89)  |  |  (ans[0]  =  =  1 2 1 ))  { 

top_flag=l;} 

else  top_flag=0; 

printf("\nDo  you  want  to  use  a  file  for  initial  conditions  (Y/N)?    "); 

scanf("%s", answer); 

samplecounter = 0; 

if((answer[0]  =  =  89)  j  |  (answer[0]  =  =  1 2 1 ))  { 

printf("\n01d  file?"); 

scanf("%s",answ); 

printf("\nEnter  name  of  file  to  be  read:    "); 

scanf("  %s", filename); 

if((fq=fopen(filename,"r"))!  =  NULL)  { 
meangamma = 0; 

fscanf(fq,"%d\n",&no_pendulums); 
fscanf(fq, "  %  1  An"  ,&alpha); 
fscanf(fq,"%lf\n",&beta); 
fscanf(fq,"%lf\n",&eta); 
fscanf(fq,"%lf\n",&omega); 
DOFOR(i,no_pendulums)  { 
fscanf(fq,"%lf   %lf  %lf  %lf\n", 
&omegaO[i],&gamma[i],&coordinate[i],&momentum[i]); 
meangamma = meangamma + gamma  [i] ; 

} 

meangamma = mean_gamma/no_pendulums ; 

fclose(fq); 

} 

else  printf(" Can't  open  file  requested."); 

}  /*  if((ans...  */ 

else  { 

printf("\nEnter  number  of  pendulums  to  use:    "); 

scanf("  %d"  ,&no_pendulums); 

printf("\nEnter  mode  amplitude:    "); 

scanf("  %  If '  ,&mode_amp); 

printf("\nEnter  modulation  amplitude:    "); 

scanf("  %  If '  ,&max_amp); 

printf("\nEnter  nonlinear  coefficient  alpha  (+/-  1  ONLY):    "); 

scanf("%lf,&alpha); 

DOFOR(k,no_pendulums)  { 

if(alpha<0)  coordinate^] =pow((-l),k)*(mode_amp 

+  max_amp*sin(2*PI*k/no_pendulums)); 

else  coordinate^]  =  mode_amp+max_amp*sin(2*PI*k/no_pendulums); 

momentum[k]=0; 

pinnedelementsfk] =0; 

}    /*  DOFOR  */ 

printf("\nEnter  coupling  coefficient  gamma:    "); 
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scanf("  %  If '  ,&mean_gamma); 
DOFOR(i,no_pendulums)  { 
gamma[i]  =  meangamma; 
omegaO[i]  =  l; 

} 

printf("\nEnter  drive  amplitude  eta:    "); 

scanf("%lf',&eta); 

printf("\nEnter  drive  frequency  omega:    "); 

scanf("  %lf '  ,&omega); 

printf("\nEnter  dissipation  constant  beta:    "); 

scanf("%lf",&beta); 

}   /*  else  */ 

printf("\nEnter  number  of  steps  per  response  cycle:  "); 

scanf("  %d"  ,&step_size); 

if  ((step_size/2)  !=  0)g  =  .5; 

else  g=0; 

time_int=2*PI/(step_size*omega); 

DOFOR(i,15)flags[i]=0; 

DOFOR(i,4096)  spectrum[i]=0; 

DOFOR(i,150)  { 

oldold_coordinate[i] =0.0; 

old_coordinate[i]  =0.0; 

} 
screengraphO; 

number_clicks=0; 

mean_omegaO=  1 .0; 
peak_amp=0; 
spectrumcounter = 0; 

/*  Phase  correction  to  start  files  in  the  proper  phase*/ 
if((answ[0]  =  =  89)  |  |  (answ[0]  =  =  121))  delta =asin((omega*beta)/eta)/2; 
else  delta = (asin((omega*beta)/eta)/2)-PI; 
model_time=0.0; 
freq_inc = omegaincrement/20; 

DINC=  .02;        /*  CHANGE  THIS  TO  CHANGE  SCALE  OF  DISPLAY  */ 
PINC=2;  /*  CHANGE  THIS  TO  CHANGE  SCALE  OF  PHASE  PLOT  */ 

}  /*  USERINIT  */ 


displaygraphicsO  { 

/*  Displays  the  lattice  sites.   The  screen  updates  sequentially  through  the 

points  after  each  round  of  motion  calculations  are  complete.  */ 

int  c,i,j,k,l,m,n; 

if(no_pendulums  <  =40)  { 
first_element=0; 
DOFOR(k,no_pendulums)  { 
n=0; 
instability  element  =  =k)&&(stability_flag=  =  1))  n=  1; 
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if((pinned_elements[k]  =  =  l)&&(disp_color=  =  1))  n  =  2; 
if((pinned_elements[k]  =  =  l)&&(disp_color=  =2))  n  =  -l; 
l=old_coordinate[k]/DISPLAY_INCREMENT; 

setcolor(O); 
_setpixel((160-5*MIDDLE_ELEMENT  +  5*k),(100+l)); 
1 = coordinate[k]/DISPLA  YINCREMENT; 
_setcolor(disp_color + n); 

_setpixel((160-5*MIDDLE_ELEMENT  +  5*k),(100+l)); 
}  /*  DOFOR  */ 
}  /*  IF  */ 

else  { 
1  =  (no_pendulums-40)/2; 
first_element=l; 
DOFOR(k,40)  { 

m=old_coordinate[l  +  k]/DISPLAY_INCREMENT; 
n=0; 

if((stability_element=  =  (l  +  k))&&(stability_flag  =  =  1))  n=  1; 
if((pinned_elements[k]  =  =  l)&&(disp_color  =  =  1))  n  =  2; 
if((pinned_elements[k]  =  =  l)&&(disp_color=  =2))  n=-l; 

setcolor(O); 
_setpixel((60  +  5*k),(100+m)); 
m  =  coordinate[l  +  k]/DISPLA  YINCREMENT; 
_setcolor(disp_color+n); 
_setpixel((60  +  5*k),(100+m)); 

}  /*  DOFOR  */ 
}     /*  ELSE  */ 
}  /*  DISPLAYGRAPHICS  */ 

displaytextO  { 

/*  The  coordinate  and  other  information  below  is  a  snapshot  of  what  the 

lattice  is  like  at  the  particular  instant  in  time  when  the  key  was  pushed.*/ 

char  message[80]; 

int  c,j,k,l,m; 

_setvideomode(_DEFAULTMODE); 

printf("Time  is  :  %lf",model_time); 

printf("     System  parameters  are:  \n"); 

printf("        Gamma  %lf",mean_gamma); 

printf("        Eta  %lf',eta); 

printf("        Omega  %lf\n", omega); 

printf("        Beta  %  If, beta); 

printf("        Alpha  %lf\n",  alpha); 

printf("\nThere  are  %d  elements  in  the  system ",no_pendulums); 

printf("\n\nPress  any  key  to  continue..."); 

c=getcharO; 

while(kbhitO==0); 

printf("Element       Position       Velocity  Element       Position       Velocity\n"); 

DOFOR(j,20)  { 
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printf("    %d         %lf  %lf  %d  %lf  %lf\n", 

j,coordinate[first_element+j],momentam[first_element+j], 

(j  +  20),coordinate[first_element  4- j + 20] , 

momentumffirstelement  +  j  +  20]); 

}  /*  DOFOR  */ 

printf("\n\nPress  "G  for  graphics,  "Q  to  quit."); 

c=getchar(); 

while(kbhitO==0); 

while(kbhitO==0); 

}   /*  display_text  */ 

process_pauseO  { 

/*  Old  method  of  saving  files.  Not  sure  if  it  has  other  uses.*/ 
int  c,i,j,k,l,mode; 

FILE  *fr; 

char  filename[30],  message[80]; 
_clearscreen(_GCLEARSCREEN); 

setvideomodeCDEFAULTMODE); 
printf("Do  you  wish  to  save  this  state?   "); 
if((c  =  getche0)=  =  121)  { 
printf("\nEnter  name  of  file  to  be  written:  "); 

scanf("  %s", filename); 

if((fr=fopen(filename,"w"))!  =  NULL)  { 

fprintf(fr,"  %d\n",no_pendulums); 

fprintf(fr,"%lf\n",alpha); 

fprintf(fr,M%lf\n",beta); 

fprintf(fr,"%lAn",eta); 

fprintf(fr,"%lf\n",omega); 

DOFOR(i,no_pendulums) 

fprintf(fr,"%lf  %lf  %lf  %lf\n", 
omegaO[i],gamma[i],old_coordinate[i],momentum[i]); 

fclose(fr); 
}  /*  ifO  */ 

else  printf(" Failed  to  open  %s\n", filename); 
}  /*  if  (0)  */ 

time_int=time_int*200/period; 

printf("\nEnter  new  time  multiple  (old  multiple  is  %\f):    ",time_int); 
scanf("  %  If  ,&time_int); 
timeint = time_int*period/200; 
if(phase_flag=  =0)  _setvideomode(_MRES4COLOR); 
else  _setvideomode(_VRES16COLOR); 
}   /*  process_pause  */ 

startfiledumpO  { 

/*  Sets  up  a  file  dump  of  the  coordinates  of  the  selected  element.   Set  up 

to  record  every  time  through  for  30  times  */ 

_clearscreen(_GCLEARSCREEN); 
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_setvideomode(_DEFAULTMODE); 

printf("Enter  number  of  element  to  be  monitored:  "); 

scanf("  %d",&chosen_element); 

if(chosen_element  >  (no_pendulums- 1 )) 

printf("Out  of  range.    No  file  dump"); 

else  filedumpflag  =  1 ; 

_clearscreen(_GCLEARSCREEN); 

if(phase_flag=  =0) 

_setvideomode(_MRES4COLOR); 

else 

_setvideomode(_VRES  16COLOR); 

counter  1=0; 

}  /*  startfiledump  */ 

stopfiledumpO  { 

/*  Dumps  the  collected  data  to  a  file.  */ 

char  filename[30]; 

int  c; 

FILE  *fs; 

_clearscreen(_GCLEARSCREEN); 

_setvideomode(wDEFAULTMODE); 

printf("\nEnter  name  of  file  to  be  written:  "); 

scanf("  %s",  filename); 

if((fs  =  fopen(filename,"w"))!  =  NULL)  { 

DOFOR(c,8000)  fprintf(fs,"%lf  %df\n",time_series_disp[c],gst[c]); 

_clearscreen(_GCLEARSCREEN); 

counter  1=0; 

if(phase_flag==0) 

_setvideomode(_MRES4COLOR); 

else  _setvideomodeCVRES16COLOR); 

}  /*  if(0)  */ 

}  /*  stopfiledump  */ 

monitorstabilityO  { 

/*  Selects  the  lattice  site  to  be  monitored  */ 

if(stability_flag==0)  { 

stability_flag=l; 

disp_color=2; 

setvideomodeCDEFAULTMODE); 
printf("Enter  number  of  element  to  be  monitored:  "); 
scanf("  %d",&stability_element); 
stabilityrestart    0; 

} 
else  { 

stab  i  1  ity  _fl  ag = 0 ; 

disp_color=l; 

} 
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if(phase_flag==0) 

_setvideomodeCMRES4COLOR); 

else 

_setvideomodeCVRES  16C0L0R); 

} 

stabilitycheckO  { 

/*  The  criterion  is  set  for  1  % .  Appears  to  check  to  see  if  20  consecutive 

peaks  are  within  the  criterion  of  each  other.  */ 

int  i,k; 

k=0; 

DOFOR(i,19)  { 

stable_fiag=  1; 

if((peak_record[i+ 1]  >  (peak_record[i]*(l  +  STABILITYINCREMENT)))  j 

(peak_record[i+l]<(peak_record[i]*(l  -  STABILITYJNCREMENT  ))))   { 

stable_flag=0; 

break; 

}  /*  if  */ 

}     /*  DOFOR  */ 

if(stable_flag=  =  l)  { 

disp_color=l; 

} 

} 

pinelementsO  { 

/*  Pins  an  element  at  rest  */ 

int  ans,i, check; 

_setvideomode(_DEFAULTMODE); 

check =0; 

printf("\nEnter  number  of  element  to  be  pinned:    "); 

scanf("%d",&ans); 

if((ans  <  0)  |  (ans  >  (no_pendulums-l)))  { 

check=l; 

} 
if(check==0)  { 

if(pinned_elements[ans]  =  =0)  { 

pinned_elements[ans]  =  1 ; 
coordinate[ans] = 0; 
momentum[ans] =0; 

} 

else  pinned_elements[ans]=0; 

}  /*  if  */ 

if(phase_flag==0) 

_setvideomode(_MRES4COLOR); 

else 

_setvideomode(_VRES  16COLOR); 

I 
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kickelementsO  { 

/*  Equivalent  to  banging  an  element  with  a  hammer.    No  control  over  when 

the  hit  occurs  in  the  cycle.  */ 

int  ans, check; 

double  amount; 

_setvideomode(_DEFAULTMODE); 

check =0; 

printf("\nEnter  number  of  element  to  kick:    "); 

scanf("%d",&ans); 

if((ans  <  0)  |  (ans  >  (no_pendulums-l)))  { 

check=l; 

} 

if(check==0)  { 

printf("\n    Enter  amount  to  kick:    "); 

scanf("  %  XT  ,&amount); 

coordinate[ans]  =  coord  inate[ans]  +  amount; 

}  /*  if  */ 

if(phase_flag==0) 

_setvideomode(_MRES4COLOR); 

else 

_setvideomode(_VRES  16COLOR); 

} 

stability_restartO  { 

/*  Changes  the  color  to  show  restart  and  zeros  the  stability  check  matrix  */ 
int  i; 

char  buffi  [50]; 
disp_color=2; 
stable_flag=0; 

DOFOR(i,20)  peak_record[i]=0; 
peak_record[  15]  =  10; 
peak_amp=0; 
sprintf (buffi, "Eta  %.31f  Omega  %.31f  Gamma  %.21f",eta,omega,mean_gamma); 

_moveto(10,10); 

_setfont("t'tms  rmn'bn5"); 

_outgtext(buffl); 

r 

save_stateO  { 

/*  Automatic  state  saver.   It's  supposed  to  keep  the  program  from  crashing 
when  a  lattice  flys  off  numerically.  */ 

int  i; 
char  filename[]  =  "savedsta"; 

FILE  *fl; 
/*  change  the  number  below  to  change  strokes  between  autosaves  */ 
if(number_clicks  =  =  100)  { 
if(chdir("E:\\MODEL\\CLEON")  !=  0)  { 
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screentextO; 

printfffailed  to  change  to  E:\\MODEL\\CLEON.\n"); 

stability  restartO;  } 
p_flag=0; 
while(p_flag==0)  { 
eqnmotionO; 

} 
if((fl  =  fopen(filename,"w"))  !=  NULL)  { 

fprintf(fl,"  %d\n",no_pendulums); 
fprintf(fl,"%lf\n",  alpha); 
fprintf(fl,"%mn",beta); 
fprintf(fl,"%lf\nH,eta); 
fprintf(fl,"%l  An", omega); 
DOFOR(i,no_pendulums) 
fprintf(fl,"%lf  %lf  %lf  %lf\n", 
omegaO[i],gamma[i],old_coordinate[i],old_momentum[i]); 
fclose(fl); 

number_clicks  =  0;  } 
else   {  screentextO; 

printf(" Failed  to  open  %s\n",filename); 
stabilityrestart;  }  } 
if(chdir("E:\\MODEL")  !=  0)  { 
screen_textO; 

printf("failed  to  change  back  to  E:\\MODEL.\n"); 
stabilityrestartO;  } 
else  + +number_clicks; 
screengraphO;  } 
/*  save_state  */ 

savestateuO  { 

/*  User  initiated  state  save.   Normal  way  to  save  lattice  data  files.   Saves 
data  for  the  lattice  at  a  turning  point  (within  one  time  step).  */ 

int  i; 

char  filename[30]; 

FILE  *fp; 

screentextO; 

printf("\nEnter  name  of  file  to  save  state  in:  "); 

scanf("  %s" , filename); 
p_flag=0; 
while(p_flag==0)  { 
eqnmotionO; 

} 

printf("Vn%lf  %lf',model_time,delta); 
if((fp=fopen(filename,"w"))  !=  NULL)  { 
fprintf(fp,"%d\n",no_pendulums); 

rprintf(fp,"%lf\n",  alpha); 

fprintf(fp,"%lf\n",beta); 
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fprintf(fp,"%lf\n",eta); 
fprintf(fp,"  %lf\n", omega); 

DOFOR(i,no_pendulums) 

fprintf(fp,"%lf   %lf  %lf   %lf\n", 
omegaO[i],gamma[i],old_coordinate[i],old_momentum[i]); 

if(top_flag=  =  l){ 

fprintf(fp,"%lf  %lf  %lf  %lAn",oldold_coordinate[stability_element], 
old_coordinate[stability_element],coordinate[stability_element],time_int); 

fprintf(fp,"%lf  %lf  %lf,,oldold_momentum[stability_element], 
old_momentum[stability_element],momentum[stability_element]);} 

fclose(fp);  } 

else  printf(" Failed  to  open  %s\n", filename); 
pause_flag=0; 

printf("\nEnter  number  of  steps  per  response  cycle:  "); 
scanf("  %d"  ,&step_size); 
if  ((step_size%2)  !=  0)  g=.5; 
else  g=0; 

time_int=2*PI/(step_size*omega); 
screengraphO;  } 
/*   save_stateu   */ 

retrievestateO  { 

/*  Retrieves  the  auto  saved  file  when  the  lattice  crashes.  */ 

int  i,c; 

char  filename[]  =  "savedsta"; 

FILE  *fg; 

screentextO; 

printf("\nLATTICE  CRASH! !!!\n"); 

printf("Do  you  want  to  retrieve  the  latest  stored  state(Y/N).\n"); 
printf("    Eta    %lf     Omega     %lAn",eta,  omega); 
if((c=getche0)=  =  121)  { 
if(chdir("E:\\MODEL\\CLEON")  !=  0)  { 
printf("Failed  to  change  directories. \n"); 

} 
} 
else  {  printf("\nEnter  name  of  file  to  retrieve.:  "); 

scanf("%s", filename);  } 
if((fg=fopen(filename,"r"))!  =  NULL)  { 

mean_gamma=0; 
fscanf(fg,"  %d\n",&no_pendulums); 
fscanf(fg,"  %lf\n",&alpha); 
fscanf(fg,"%lf\n",&beta); 
fscanf(fg,"%lAn",&eta); 
fscanf(fg, "  %lAn"  ,&omega); 
DOFOR(i,no_pendulums)  { 
fscanf(fg,"%lf  %lf  %lf  %lf\n", 
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&omegaO[i],&gamma[i],&coordinate[i],&momentum[i]); 
meangamma = meangamma + gamma[i] ; 

} 

fclose(fg); 

counter2=0; 

delta=atan((omega*beta)/sqrt(SQR(eta)-(SQR(omega)*SQR(beta))))/2; 

model_time=0.0; 

} 

else  printf("Can't  open  file  requested."); 

timeint = time_int*200/period ; 

printf("\nEnter  new  time  multiple  (old  multiple  is  %lf):    ",time_int); 

scanf( "  %  If  ,&time_int) ; 

time_int=time_int*period/200; 

ifCchdirC'EiWMODEL")  !=  0)  { 

screentextO; 

printf("failed  to  change  back  to  E:\\MODEL.\n");  } 

screengraphO; 

} 
/*  retrieve_state  */ 

screentextO  { 

_clearscreen(_GCLEARSCREEN); 
textflag  =  1 ; 
graphics_flag=0; 
spectrum  flag =0; 

_setvideomode(_DEFAULTMODE);  } 
/*  screentext   */ 

screengraphO  { 

_clearscreen(_GCLEARSCREEN); 
_setvideomode(wMRES4COLOR) ; 

setcolor(l); 
phase_flag=0; 
spectrumflag = 0; 
text_flag=0; 
graphics_flag=l;  } 
/*  screen_graph    */ 

linearkickO  { 

/*  Lattice  perturbation  up  to  a  maximum  of  3%.   Sort  of  random  as  to  sign  of 
kick  direction  for  each  site.   Uses  the  %  of  the  maximum  amplitude  site.  */ 
int  aj,i; 
double  k,b; 

struct  dostimet  time2; 
_dos_gettime(&time2); 
srand(time2.hsecond); 
DOFOR(i,no_pendulums)  { 
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b  =  rand(); 

/*  change  33.33  to  change  %,  now  3%  */ 

k=b/(33.33  *  32767); 

ampkick = k*peak_amp; 

if(b  >  16384)  { 

coordinate[i]  =  coordinate^]  +  amp_kick; 

} 

else  coordinated  =  coordinate[i]  -  ampkick; 

} 

peak_amp=0.0; 
}  /*  linearkick  */ 

timewaitO  { 

/*  part  of  a  routine  to  get  the  lattice  at  the  turning  point.  */ 
t_flag=0; 

while(t_flag==0){ 

eqnmotionO; 

} 

}/*  timewait  */ 

init_phasplotO  { 

/*  Sets  up  the  phase  plot  to  look  at  the  phase  plane  of  sites.  */ 

int  c,i,j,k,l; 

char  ans[5],message[40],txt[3]; 

float  dummy; 

spectrum_flag=0; 

_setvideomode(_DEFAULTMODE); 

DOFOR(i,5)  phase_elements[i]=0; 

printf("\nDo  you  want  Poincare  sections?   "); 

if((c=getchO)=  =  121)  phase_flag=2; 

else  phase_flag=l; 

printf("\n Which  elements  do  you  wish  to  monitor  (999  to  finish):  "); 

i=-l; 

k=0; 

while((i+  +  !=999)&&(k<5))  { 

scanf("%d",&i); 

phase_elements[k  +  +  ]  =  i  + 1 ; 

}  /*  while  */ 

DOFOR(i,5)  { 

if(phase_elements[i]=  =  1000)  { 

phase_elements[i] =0; 

} 

} 

_setvideomode(_VRES  16COLOR); 

_clearscreen(_GCLEARSCREEN); 

graphics_flag=0; 
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_settextcolor(7); 
colors[0]=2; 
colors[l]  =  3; 
colors[2]=9; 
colors[3]  =  14; 
colors[4]  =  13; 
DOFOR(i,120)  { 
_setpixel(318,4*i); 

} 

DOFOR(i,160)  { 
_setpixel(4*i,238); 

} 

D0F0R(i,6)  { 

_setpixel(3 19, 10*SCREEN_CORRECTION  FACTOR + 1 

+  (80/SCREEN_CORRECTION_FACTOR)*(i+l)); 

} 

DOFOR(i,8)  { 

_setpixel(80*(i+l),239); 

} 

_moveto(480,450); 
DOFOR(i,5)  { 
if((phase_elements[i])!=0)  { 
c = phase_elements[i]- 1 ; 

/*         itoa(c,txt,10); 

outtext(txt);  */ 
printf("%d   ",c); 
_setcolor(colors[i]); 
DOFOR(j,4)  { 
DOFOR(k,8)  { 

_setpixel(550+  10*i  +  k,10+j); 
} 
} 
} 

}   /*  DOFOR  */ 
printf("       LAT2X486.C"); 

printf("\nGamma:    %lf  Eta:    %lf  Beta:    %lf",mean_gamma,eta,beta); 
printf("\nOmega:    %lf  Alpha:    %lf    Time  interval:    %lf, 
omega,alpha,time_int); 
dummy  =  l/(time_int*omega); 
phaseint = dummy/ 1 ; 
phase_counter=0; 
} 

phasplotO  { 

/*  Displays  the  phase  for  chosen  elements.  */ 

int  i,j,k,l; 

double  speed,  position,fast_coordinate,fast_momentum; 
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if(phase_flag  =  =  1 )  { 

D0F0R(i,5)  { 

if((j=phase_elements[i]-l)!=-l)  { 

speed  =  PHASE_INCREMENT*momentum[j]; 

_setcolor(colors[i]); 

position =PHASE_INCREMENT*coordinate[j]; 

fast_coordinate=(position*cos(omega*model_time) 

-speed*sin(omega*model_time)/omega) 

/(3*SCREEN_CORRECTION_FACTOR); 

fastmomentum  =  (position*sin(omega*model_time) 

+  speed  *cos(omega*model_time)/omega)/4; 

_setpixel((320+320*fast_coordinate),(240+240*fast_momentum)); 

}  /*  if(0)  */ 

}     /*  DOFOR  */ 

}       /*  if  */ 

else  if(phase_flag=  =2)  { 

if(phase_counter  =  =  phaseint)  { 
phase_counter=0; 
DOFOR(i,5)  { 

if((j=phase_elements[i]-l)!  =  -l)  { 
speed  =  PH  ASE_INCREMENT*momentum[j] ; 
_setcolor(colors[i]); 

position  =  PH  ASE_INCREMENT*coordinate[j] ; 
fastcoordinate = (position*cos(omega*model_time) 
-speed*sin(omega*model_time)/omega); 
fastmomentum  =  (position*sin(omega*model_time) 
+  speed*cos(omega*model_time)/omega)/4; 
_setpixel  ((320 + 320  *fast_coord  inate) ,  (240 + 240  *f fastmomentum)) ; 
}  /*  if(0)  */ 
}     /*  DOFOR  */ 

}  /*  if  */ 
}  /*  else  ifO  */ 
}  /*  phasplot  */ 

/*  FFT  Routines  for  lattice  programs  */ 

computespectrumOl  /*  Uses  algorithm  from  p.  411  of  Numerical  Recipes  in  C  */ 

int  n,mmax,m,j,istep,i,nn,temp; 

double  wtemp,wr,wpr,wpi,wi,theta,tempr,tempi; 

nn  =  2048; 

n  =  nn<  <  1; 

DOFOR(i,2048)  { 
nn=(i«&  1024)/ 1024; 
nn  =  nn  +  (i&512)/256; 
nn=nn+(i&256)/64; 
nn=nn+(i&128)/16; 
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nn=nn+(i&64)/4; 

nn=nn+(i&32); 

nn=nn+(i&16)*4; 

nn=nn+(i&8)*16; 

nn  =  nn  +  (i&4)*64; 

nn=nn+(i&2)*256; 

nn=nn+(i&l)*1024; 

temp2  [2  *nn]  =  spectrum[2  *i] ; 

temp2[2*nn+  l]  =  spectrum[2*i+ 1]; 

} 

DOFOR(i,4096)  spectrum[i]=temp2[i]; 

mmax=2; 

while(n>mmax)  { 

istep  =  2*mmax; 

theta=2*PI/mmax; 

wtemp  =  sin(0.5*theta); 

wpr = -2 . 0  *  wtemp  *  wtemp ; 

wpi=sin(theta); 

wr=1.0; 

wi=0.0; 

for(m=0;m<(mmax-l);m+=2)  { 

for(i  =  m;i  <  =  n;i  +  =  istep)  { 

j  =  i  +  mmax; 

tempr=wr*spectrum[j]-wi*spectrum[j  +  l]; 

tempi  =  wr*spectrum(j  + 1  ]  +  wi*spectrum[j]; 

spectrum[j] =spectrum[i]-tempr; 

spectrum[j+  l]  =  spectrum[i  +  l]-tempi; 

spectrum[i]+  =tempr; 

spectrum[i+ 1]+  =  tempi; 

} 

wr = (wtemp = wr)*wpr-wi*wpi + wr ; 

wi  =  wi*wpr + wtemp  *wpi  +  wi; 

} 

mmax  =  istep; 

} 
} 

plotspectrumO  { 

int  i  j; 

double  freq,magnitude,mag,max; 

int  c; 

FILE  *fp; 

if(dump_spectrum_flag=  =  l)  { 

printf("\nDumping  spectrum  now"); 

fp = fopen( "  spectrum .  out " , "  w " ) ; 

for(i=l;i<2048;i++)  { 

if(time_int!  =0)  freq=i/(2048*time_int); 
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else  printf("Time_int  was  zero!"); 
magnitude=sqrt(spectrum[(2*i)]*spectrum[(2*i)] 
+spectrum[(2*i+  l)]*spectrum[(2*i+ 1)]); 
fprintf(fp,"%lf  %lf  \n",freq,magnitude); 

} 

fclose(fp); 

dumpspectrumflag = 0; 

} 
_clearscreen(J3CLEARSCREEN); 

max=0; 

for(i=l;i<2049;i++)  {  /*  this  loop  calculates  the  max  spectral  component*/ 

magnitude=sqrt(spectrum[(2*(i))]*spectrum[(2*(i))] 

+spectrum[(2*(i)+  l)]*spectrum[(2*(i)  + 1)]); 

if(magnitude  >  max)  max  =  magnitude; 

} 

_setcolor(3); 

DOFOR(i,640)  _setpixel(i,460); 
DOFOR(i,l  1)  _setpixel(51*i,461),_setpixel(51*i,462); 
printf("Spectrum  for  element  %d       LATTIC2X.C",spectrum_element); 
printf("\nAlpha:  %lf  Beta:  %lf  Gamma:  %lf  Eta:  %lf  Omega:  %lf\ 
alpha,beta,mean_gamma,eta,omega); 
printf("\nMax  amp  =  %lf       ",max); 
freq = 5 1 2/(2048*time_int) ; 

printf("Max  freq  shown:  %lf      Time  interval:    %lf',freq,time_int); 
setcolor(l); 
DOFOR(i,512)  { 
j=0; 

magnitude=sqrt(spectrum[2*i]*spectrum[2*i] 
+  spectrum[2*i+  l]*spectrum[2*i+ 1]); 
mag = magnitude  *400/max; 
while(j+  +  <mag)  { 
_setpixel(i,460-j); 
}  /*  while  */ 
}    /*  DOFOR  */ 
} 

initwaterfallO  { 

/*  Sets  up  waterfall  type  display  to  see  trends  in  latice  motion.  Good  for 
observing  kink  or  breather  drift.  */ 
int  c,i,j,k,l,m; 
char  ans[5]; 

_setvideomode(_VRES  16COLOR); 
setcolor(l); 
DOFOR(i,600)  _setpixel(20+i,20); 
DOFOR(i,450)  _setpixel(20,i+20); 
DOFOR(i,61)  { 
if((i=  =  10)|  (i=  =20)|(i=  =32)|  (i=  =42)  |  (i=  =52))_setcolor(2); 
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D0F0R(j,113)  { 
_setpixel(20+  10*i,20+4*j); 

} 

setcolor(l); 

} 

_setcolor(2); 

DOFOR(i,ll)  { 
DOFOR(j,160)  { 
_setpixel(20+4*j,16+40*(i  + 1)); 

} 

} 

waterfallcounter = 0; 

color_counter=2; 
} 

dispwaterfallO  { 
int  c,i  j,k,l; 

if(waterfall_counter+  +  <  1 10)  { 
if(color_counter+  +  =  =  14)  { 
color_counter=2; 
setcolor(colorcounter); 

} 

else  { 
setcolor(colorcounter); 

} 

DOFOR(i,no_pendulums)  { 
if(waterfall_counter<55)  { 
_setpixel(20+2*i,24+8*waterfall_counter-7*wat_inc*(coordinate[i])); 

} 

else  { 
_setpixel(340+2*i,24+8*(waterfall_counter-55)-7*wat_inc*(coordinate[i])); 

} 

} 

} 

else  { 

wait_flag=l; 

} 

} 

eqnmotionO  { 

/*     The  following  loop  calculates  a  round  of  velocities  for  one 
iteration  of  the  model  clock.   This  is  where  the  equations 
of  motion  need  to  be  located,  if  one  wishes  to  change  the 
model's  physics!    */ 
int  i,j,k; 

DOFOR(i,no_pendulums) 
if(coordinate[i]  >  100)  { 
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retrieve_state(); 

} 
DOFOR(i,no_pendulums)  { 

if(i==0){ 

acceleration = gamma[i]  *(coordinate[i  + 1  ]  4-  coordinate[no_pendulums- 1  ] 

-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])+2*eta*cos(2*omega*model_time-delta))*coordinate[i] 

+  alpha*CUB(coordinate[i]); 

} 

else  if(i==(no_pendulums-l))  { 

acceleration=gamma[i]*(coordinate[i-l]  +  coordinate[0] 

-2*coordinate[i])-beta*mornentum[i] 

-(SQR(omegaO[i])  +  2*eta*cos(2*omega*model_time-delta))*coordinate[i] 

+  alpha*CUB(coordinate[i]); 

} 

else  { 

acceleration =gamma[i]*(coordinate[i+  l]  +  coordinate[i-l] 

-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])  +  2*eta*cos(2*omega*model_time-delta))*coordinate[i] 

+alpha*CUB(coordinate[i]); 

} 

oldold_momentum[i] = old_momentum[i] ; 

old_momentum[i]  =  momentum[i]; 

momentum[i]  =  momentumfi]  +  acceleration*time_int; 

}   /*  DOFOR  */ 

DOFOR(i,no_pendulums)  { 

oldold_coordinate[i]=old_coordinate[i]; 

old_coordinate[i]  =  coordinate[i]; 

if(pinned_elements[i]==0)  { 

coordinate^]  =  coordinate[i]  +  momentum[i]  *time_int; 

} 

/*  Records  data  for  file  dump  */ 

if((i=  =chosen_element)&&(file_dump_flag=  =  1))  { 

/*   if(counter2+  +  =  =30)  {  */ 
time_series_disp[counterl  +  +]=coordinate[i]; 
gst[counterH- +]=g; 

/*       counter2=0; 

}*/ 

if(counterl  =  =  8000)  stopfiledumpO; 

}  /*  if((i...  */ 

/*  Stores  energy  information  */ 

if(energy_flag=  =  1)  { 

energy  =  energy + 0.5  *(SQR(momentum[i]) 

+  gamma[i]*SQR((coordinate[i+  l]-coordinate[i])) 

+  SQR(coordinate[i]))-alpha/4*pow(coordinate[i],4); 
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}  /*  if(energy...)  */ 

/*  Stores  FFT  information  */ 
if  ((spectrum  flag  =  =  l)&&(spectrum_element=  =i) 
&&((sample_counter+  +)=  =  1))  { 
if(spectrum_counter  <  2048)  { 
spectrum[2*spectrum_counter]  =  coordinate[i]; 
spectrum[2*spectrum_counter+ 1]=0; 
spectrum_counter+  + ; 
samplecounter = 0; 

} 

else  { 

computespectrumO ; 

_clearscreen(_GCLEARSCREEN); 

_setvideomode(_VRES  16COLOR); 

plotspectrumO; 

phase_flag=0; 

text_flag=0; 

graphics_flag=0; 

DOFOR(j,4096)  spectrum[j]=0; 

spectrumcounter = 0; 

samplecounter = 0; 

}  /*  else  */ 

}  /*  if  ((spec...  */ 

/*  Records  the  peaks  for  the  stability  checker  */ 

if(((i=  =stability_element)&&(stability_flag=  =  l)&&(stable_flag=  =0))  j 

((i=  =(no_pendulums-l))&&(waterfall_flag=  =  1)))  { 

if((coordinate[i]  >  0)&&(coordinate[i]  <  old_coordinate[i]) 

&&(peak_flag==0))  { 

peak_flag=l; 

if(pause_flag=  =  1)  pause_flag2=  1; 

DOFOR(k,19)  peak_record[k]=peak_record[k+ 1]; 

peakrecord  [19]  =  coord  inate[i] ; 

if(stability_flag=  =  1)  stabilitycheckO; 

if(waterfall_flag=  =  1)  dispwaterfallO; 

} 

else  if(coordinate[i]  <  0)  peak_flag=0; 

} 

/*  Records  peak  amplitude  for  display  and  linear  kick  */ 
if((coordinate[i]  >  0)&<&(coordinate[i]  <  oldcoordinatefi]))  { 
if(old_coordinate[i]  >  peak_amp){ 
peakamp  =  old_coordinate[i]; 

} 
} 
if((i=  =stability_element)&&(stability_flag=  =  l)&&(stable_flag=  =  1)){ 
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if(((coordinate[i]-old_coordinate[i])  >  =  0)  && 

((old_coordinate[i]-oldold_coordinate[i])  <  =  0)  &&  (coord inate[i]  <  0)){ 
p_flag=l; 

} 
} 

}   /*  DOFOR  */ 

/*  Increments  time  and  resets  it  at  2Pi.   Causes  problems  in  eqns  of 

motion  if  not  reset  */ 

modeltime = modeltime + timeint; 

if((g+l)  ==  (step_size/2)){ 

model_time=0.0; 

t_flag=l; 

g=0; 

} 

else  g+  +  ; 

/*  Displays  energy  information  gathered  above.  */ 

if(graphics_flag=  =  1)  display_graphics(); 

if((text_flag=  =  l)&&(energy_flag=  =  1))  { 

printf("\n  Total  Energy  at  time  %lf  is  %lf',model_time,energy); 

energy_counter+  + ; 

if(energy_counter==20)  { 

printf("\n\nStrike  any  key  to  continue..."); 

while(kbhitO==0); 

energycounter =0; 

} 
} 
if(phase_flag!=0)  { 

phasplotO; 

phase_counter+  +  ;       /*  Used  for  Poincare  sections  */ 

} 
} 

2.        Upper  Cutoff,  Global  Drive,  Equation  of  Motion 

eqnmotionO  { 
/*      The  following  loop  calculates  a  round  of  velocities  for  one 
iteration  of  the  model  clock.   This  is  where  the  equations 
of  motion  need  to  be  located,  if  one  wishes  to  change  the 
model's  physics!   */ 
int  i,j,k; 

DOFOR(i,no_pendulums) 
if(coordinate[i]  >  100)  { 
retrievestateQ; 
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} 

DOFOR(i,no_pendulums)  { 
if(i==0){ 

acceleration=gamma[i]*(coordinate[i+l]  +  coordinate[no_pendulums-l] 
-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])  +  2*eta*cos(2*omega*model_time-delta))*coordinate[i] 
+  alpha*CUB(coordinate[i]); 

} 

else  if(i==(no_pendulums-l))  { 

acceleration =gamma[i]  *(coordinate[i- 1  ]  +  coord  inate[0] 

-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])+2*eta*cos(2*omega*model_time-delta))*coordinate[i] 

+  alpha*CUB(coordinate[i]); 

} 

else  { 

acceleration = gammafi]  *(coordinate[i  + 1  ]  +  coordinate^- 1  ] 

-2*coordinate[i])-beta*momentum[i] 

-(SQR(omegaO[i])  +  2*eta*cos(2*omega*model_time-delta))*coordinate[i] 

+ alpha*CUB(coordinate[i]); 

} 

oldmomentumfi]  =  momentumf  i]; 

momentum[i]  =  momentum[i]  +  acceleration*time_int; 

}   /*  DOFOR  */ 

DOFOR(i,no_pendulums)  { 

oldoldcoordinatefi] = old_coordinate[i] ; 

old_coordinate[i]  =  coordinate[i]; 

if(pinned_elements[i]==0)  { 

coordinate[i]  =  coordinate^]  -I-  momentum[i]  *time_int; 

} 

3.        Lower  Cutoff,  End  Driven,  Equation  of  Motion 

eqnmotionO  { 
/*      The  following  loop  calculates  a  round  of  velocities  for  one 
iteration  of  the  model  clock.   This  is  where  the  equations 
of  motion  need  to  be  located,  if  one  wishes  to  change  the 
model's  physics!    */ 
int  i,j,k; 

DOFOR(i,no_pendulums) 
if(coordinate[i]  >  1000)  { 
retrievestateO; 

} 

DOFOR(i,no_pendulums)  { 
if(i==0){ 

acceleration =gamma[0]*(coordinate[l]-coordinate[0]) 
-beta*momentum[0] 
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-((coordinate[i]/sqrt(l  +  2*SQR(coordinate[i])))  +  eta*cos(omega*model_time)); 

} 
else  if(i==no_pendulums-l){ 

acceleration = gamma[i]  *(coordinate[i- 1  ]-coordinate[i]) 

-beta*momentum[  i] 

-(coordinate[i]/sqrt(l  +  2*SQR(coordinate[i]))); 

} 
else  { 

acceleration = gamma[i]  *(coordinate[i  + 1  ]  +  coordinate^- 1  ] 
-2*coordinate[i])-beta*momentum[i] 
-(coordinate[i]/sqrt(l  +  2*SQR(coordinate[i]))); 

} 
old_momentum[i]  =  momentum[i]; 

momentum[i]  =  momentum[i]  -I-  acceleration*time_int; 
}   /*  DOFOR  */ 

DOFOR(i,no_pendulums)  { 
old_coordinate[i]  =  coordinate[i]; 
if(pinned_elements[i]=  =0)  { 
coordinate[i]  =  coord  inate[i]  +  momentum[i]*time_int; 

} 

4.        Upper  Cutoff,  End  Driven,  Equation  of  Motion 

eqnmotionO  { 
/*      The  following  loop  calculates  a  round  of  velocities  for  one 
iteration  of  the  model  clock.   This  is  where  the  equations 
of  motion  need  to  be  located,  if  one  wishes  to  change  the 
model's  physics!    */ 
int  i,j,k; 

DOFOR(i,no_pendulums) 
if(coordinate[i]  >  100)  { 
retrievestateO; 

} 

DOFOR(i,no_pendulums)  { 
if(i==0){ 

acceleration = gamma[0]  *(coordinate[  1  ]-coordinate[0]) 
-beta*momentum[0] 
-(SQR(omegaO[0])*coordinate[0]  +  eta*cos(omega*model_time))  +  alpha*CUB(coordinate[0]); 

} 

else  if(i==no_pendulums-l){ 

acceleration =gamma[i]*(coordinate[i-l]-coordinate[i]) 

-beta*momentum[i] 

-(SQR(omegaO[i])*coordinate[iJ)  +  alpha*CUB(coordinate[i]); 

} 

else  { 

acceleration = gamma[i]  *(coordinate[i  + 1  ]  +  coordinate[i- 1  ] 
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-2*coordinate[i])-beta*momentum[i] 
-(SQR(omegaO[i])*coordinate[i])  +  alpha*CUB(coordinate[i]); 

} 

old_momentum[i] =momentum[i]; 

momentum[i]  =  momentum[i]  +  acceleration  *time_int; 

}   /*  DOFOR  */ 

DOFOR(i,no_pendulums)  { 

old_coordinate[i]= coordinate^]; 

if(pinned_elements[i]==0)  { 

coordinate^]  =  coordinate[i] + momentum[i]  *time_int; 

} 
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